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ABSTRACT 

A  fractional-step  splitting  scheme  breaks  the  full  Navier-Stokes  equations 
into  explicit  and  implicit  portions  amenable  to  the  calculus  of  variations. 
Beginning  with  the  functional  forms  of  the  Poisson  and  Helmholtz  equations,  we 
substitute  finite  expansion  series  for  the  dependent  variables  and  derive  the 
matrix  equations  for  the  unknown  expansion  coefficients.  This  method  employs  a 
new  splitting  scheme  which  differs  from  conventional  three-step  (non-linear, 
pressure,  viscous)  schemes.  The  non-linear  step  appears  in  the  conventional, 
explicit  manner,  the  difference  occurs  in  the  pressure  step.  Instead  of  solving  for 
the  pressure  gradient  using  the  non-linear  velocity,  we  add  the  viscous  portion  of 
the  Navier-Stokes  equation  from  the  previous  time  step  to  the  velocity  before 
solving  for  the  pressure  gradient.  By  combining  this  "predicted"  pressure 
gradient  with  the  non-linear  velocity  in  an  explicit  term,  and  the  Crank- 
Nicholson  method  for  the  viscous  terms,  we  develop  a  Helmholtz  equation  for  the 
final  velocity. 


LIST  OF  SYMBOLS 

a  =  value  of  non-homogeneous  essential  boundary  condition 

g  as  value  of  non-homogeneous  natural  boundary  condition 

A  =  cross-sectional  area 

dA  =  differential  surface  area  vector 

Bo  =  23/12,  coefficient  used  in  Adams-Bashforth  method 
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Bt  =-16/12,  coefficient  used  in  Adams-Bashforth  method 
B2  =  5/12,  coefficient  used  in  Adams-Bashforth  method 
Bi£c  —  coefficients 
(Cij)ibc  =  coefficients 

Dn  =  Legendre  collocation  derivative  of  order  n 

Dij  =  derivative  of  expansion  polynomial  j  evaluated  at  node  i 

ex,  ey,  ez  =  unit  vectors  in  the  direction  of  the  co-ordinate  axes 

f  =  body  force  vector 

E  =  total  number  of  elements 

M  =  moment  vector 

F  =  — A2Vn*1,  also  represents  a  force  vector 
J[P]  =  a  functional  depending  on  P 
L  =  periodic  length  of  domain 
L2  =  square-integrable  function 
Lk(x)  =  Legendre  polynomial 

it,  jt,  kt  =  the  maximum  number  of  nodes  in  the  r,  s,  t  directions 
respectively 

In  =  a  system  of  interpolating  polynomials  of  order  n 

Pe  =  surface  integral  for  side  number  p  and  element  number  e 

Se  =  total  surface  integral  for  element  number  e 

[J]  =  Jacobian  matrix  for  the  transformation  between  co-ordinate  systems 
n  =  surface  unit  normal  vector 
8  =  surface  unit  tangent  vector 
R  =  position  vector 
p  =  pressure 
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pi(x)  =  infinite  sequence  of  orthogonal  functions 
p(x),  q(x),  w(x)  =  functions  in  Sturm-Liouville  equation 
P  =  p/p  +  JV-V,  dynamic  pressure 
Pn  =  a  system  of  orthogonal  polynomials  of  degree  n 
S  =  an  infinite  system  of  orthogonal  polynomials 
t  =  time 

At  =  time-step  size 

t(n)  =  force  per  unit  area  on  a  surface  with  unit  normal  n 
Ui  =  discrete  expansion  coefficient  for  the  infinite  series 
iii  =  discrete  expansion  coefficient  for  the  finite  series 

V  =  velocity  vector 

V  =  velocity  after  the  non-linear  step 

V  =  velocity  after  the  explicit  viscous  step 

V  =  velocity  after  the  pressure  step 

V  =  corrected  velocity  after  the  explicit  viscous  step 
V0ut  =  average  outflow  velocity  correction 
U,  V,  W  =  x,  y,  z  velocities  respectively 
wi  =  weight  function,  integral  of  the  expansion  polynomial  over  domain 
x,  y,  z  =  coordinates  in  global  or  physical  space 

xr  =  partial  derivative  of  the  global  co-ordinate  x  with  respect  to  the  local 
co-ordinate  r 

r,  s,  t  =  coordinates  in  local  or  transformed  space 

rx  —  partial  derivative  of  the  local  co-ordinate  r  with  respect  to  the  global 
co-ordinate  x 
Greek  and  other  Symbols 
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A 


a  —  constant  in  homogeneous  boundary  condition 
0  =  constant  in  homogeneous  boundary  condition 
7  =  coefficient  in  expansion  polynomial  relationships 
6  =  differential  of  a  quantity 
$ij  =  Kronecker  delta 

^i(x)  =  a  combination  of  Legendre  polynomials 

c  =  infinitesimal  quantity 

eijk  =  alternating  tensor 

fl  =  domain  under  consideration 

dQ  —  boundary  of  domain  under  consideration 

n  =  nudn 

$  =  \p*1 
A2=  2/(i/At) 

A  =  eigenvalue  in  Sturm-Liouville  equation 

A  =  matrix  of  coefficients  representing  A  in  the  expression  Ax=b 

II  =  matrix  of  coefficients  representing  b  in  the  expression  Ax=b 

p  —  fluid  density 

a  =  stress  tensor 

a  —  element  surface  area 

°ij  =  i,  j  component  of  the  stress  tensor 

X  =  moment  arm 

V  =  element  volume 

p.  =  fluid  viscosity 

u  —  fluid  kinematic  viscosity 

<  -  Cx€x  +  Cyey  +  (zez,  vorticity  vector 
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Subscripts 

o  -  initial  value 

as  =  free-stream  value 

e  =s  element  number 

n  =  direction  normal  to  the  surface 

s  =  direction  tangential  to  the  surface 

ijk  =  co-ordinate  system  indices 

in  =  value  at  inlet 

out  =  value  at  outlet 

wall  =  value  at  wall 

x,y,z  =  streamwise,  vertical,  and  spanwise  values  respectively,  may  also 
refer  to  partial  derivatives  with  respect  to  the  global  co-ordinates  x,  y,  z 
Superscripts 

n  =  time  step  number 

e  —  element  number 

i>2, 3, 4, 5, s  =  element  side  numbers 

INTRODUCTION 

The  new  splitting  scheme,  developed  by  Wessel  (1992),  contains  variations 
on  the  original  three-step  splitting  method  proposed  by  Korczak  and  Patera 
(1986).  In  the  previous  scheme,  the  non-linear,  pressure,  and  viscous  terms  in  the 
incompressible  Navier-Stokes  equations  appear  in  separate  fractional  steps.  By 
introducing  intermediate  velocities,  solutions  of  these  equations  yield, 
consecutively,  a  velocity  field  based  on  the  non-linear,  the  pressure  and  non¬ 
linear,  and  the  viscous,  pressure,  and  non-linear  terms.  The  final  step  producing 
the  true  velocity  field. 
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The  velocity  field  resulting  from  the  non-linear  step  satisfies  no  boundary 
conditions  nor  the  incompressibility  constraint.  This  velocity  field  supplies  the 
forcing  function  for  the  Poisson  equation  for  pressure  after  applying  the 
divergence  operator  to  the  pressure  step.  The  intermediate  velocity  contained  in 
the  pressure  step  must  satisfy  the  divergence  free  constraint,  thus  it  vanishes 
from  the  Poisson  equation  for  pressure.  Instead  of  solving  a  second-order  Poisson 
equation  for  pressure,  a  first-order  equation  for  pressure  gradient  is  solved  using 
methods  from  the  calculus  of  variations-the  velocity  field  follows  directly  from 
the  pressure  step.  Inviscid  boundary  conditions  on  velocity  determine  the 
pressure  boundary  conditions;  hence,  errors  of  0(at)  occur  near  solid  boundaries. 
Finally,  the  viscous  step  employs  a  Crank-Nicholson  scheme  yielding  a  Helmholtz 
equation  with  a  forcing  function  determined  by  the  velocity  after  the  pressure 
step.  Once  again  a  variational  form  of  the  governing  second-order  differential 
equation  reduces  the  order  by  one.  The  velocity  must  satisfy  the  full,  viscous 
boundary  conditions;  however,  it  does  not  satisfy  the  incompressibility 
constraint. 

The  new  method  varies  slightly  from  the  old.  Instead  of  solving  for  the 
pressure  gradient  using  the  non-linear  velocity,  we  include  the  viscous  term  from 
the  previous  time  step.  Thus,  the  forcing  function  appearing  in  the  Poisson 
equation  for  pressure  contains  contributions  from  both  non-linear  and  viscous 
terms.  The  resulting  pressure  gradient  is  not  solved  for  the  velocity  after  the 
pressure  step;  instead,  it,  along  with  the  velocity  from  the  non-linear  step,  and  a 
Crank-Nicholson  method  for  the  viscous  terms,  produce  a  Helmholtz  equation  for 
the  full  velocity.  The  boundary  conditions  and  solution  procedure  remain 
identical  to  the  original  method.  The  boon  comes  from  including  the  viscous 
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terms  in  the  pressure  gradient  prediction,  resulting  in  a  quicker  solution  of  the 
pressure  step. 
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CHAPTER  I 


SPECTRAL  APPROXIMATION 

1.1  Spectral  Theory 

The  expansion  of  a  function  u  in  terms  of  an  infinite  sequence  of 
orthogonal  functions  {pi},  u  =  3??sjd  Uipi,  underlies  many  numerical  methods 
of  approximation.  The  most  familiar  approximation  results  apply  to  periodic 
functions  expanded  in  Fourier  series.  In  this  case,  the  i— th  coefficient  of  the 
expansion  decays  faster  than  any  inverse  power  of  i  for  smooth  functions  with 
periodic  derivatives.  The  rapid  decay  of  the  coefficients  implies  that  the  Fourier 
series  truncated  after  a  few  terms  represents  a  good  approximation  to  the 
function.  This  characteristic  refers  to  the  "spectral  accuracy"  of  the  Fourier 
method. 

Spectral  accuracy  for  smooth  but  non-periodic  functions  occurs  with  the 
proper  choice  of  expansion  functions.  Not  all  orthogonal  expansion  functions 
provide  high  accuracy;  however,  the  eigenfunctions  of  a  singular  Sturm-Liouville 
operator  allow  spectral  accuracy  in  the  expansion  of  any  smooth  function,  with 
no  a  priori  restriction  on  the  boundary  beMvinr. 

The  expansion  in  terms  of  an  orthogonal  system  introduces  a  linear 
transformation  between  u  and  the  sequence  of  its  expansion  coefficients  {iii}, 
called  the  finite  transform  of  u  between  physical  space  and  spectral  space.  Since 
the  expansion  coefficients  depend  on  all  the  values  of  u  in  physical  space,  they 
rarely  get  computed  exactly;  instead,  a  finite  number  of  approximate  expansion 
coefficients  result  from  using  the  values  of  u  at  a  finite  number  of  selected 
points— the  nodes.  This  procedure  defines  a  discrete  transform  between  the  set  of 
values  of  u  at  the  nodes  and  the  set  of  approximate,  or  discrete  coefficients. 
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With  a  proper  choice  of  nodes  and  expansion  functions,  the  finite  series  defined 
by  the  discrete  transform  represents  the  interpolation  of  u  at  the  nodes. 
Maintaining  spectral  accuracy  when  replacing  the  finite  transform  with  the 
discrete  transform  allows  use  of  the  interpolation  series  instead  of  the  truncated 
series  in  approximating  functions. 

1.2  Sturm-Liouville  Problems 

The  importance  of  Sturm-Liouville  problems  for  spectral  methods  lies  in  the 
fact  that  the  spectral  approximation  of  the  solution  of  a  differential  problem 
often  occurs  as  a  finite  expansion  of  eigenfunctions  of  a  suitable  Sturm-Liouville 
problem.  The  general  form  of  the  Sturm-Liouville  problem  satisfies 

“Hx^Pdx^  +  QU  =  Awu  in  H  e  (-1,1).  (1.2.1) 

The  real-valued  functions,  p(x),  q(x),  and  w(x),  must  behave  properly:  p(x) 
must  be  continuously  differentiable,  strictly  positive  in  (—1,1)  and  continuous  at 
x=±l;  q(x)  must  be  continuous,  non-negative  and  bounded  in  (-1,1);  the 
weight  function  w(x)  must  be  continuous,  non-negative  and  integrable  over 
(—1,1).  The  Sturm-Liouville  problems  of  interest  in  spectral,  methods  allow  the 
expansion  of  an  infinitely  smooth  function  in  terms  of  their  eigenfunctions  while 
guaranteeing  spectral  accuracy. 

1.3  Orthogonal  Systems  of  Polynomials 

Consider  the  expansion  of  a  function  in  terms  of  a  system  of  orthogonal 
polynomials  of  degree  less  than  or  equal  to  n,  denoted  by  Pn.  Assume 
{pk}k=o,i---  represents  a  system  of  algebraic  polynomials  (with  degree  of  pk=k) 
mutually  orthogonal  over  the  interval  (—1,1)  with  respect  to  a  weight  function 
w(x).  The  orthogonality  condition  requires 

J*  iPk(x)ptB(x)w(x)dx  =  0  whenever  m  ^  k.  (1.3.2) 
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The  formal  seric*.  of  a  squared ntegrable  function,  u€L2(— 1,1)$,  in  terms  of  the 
system  appears  as 

Su  =  £  UkPk(x),  (1.3.3) 

k  =  0 

when  the  expansion  coefficients  Uk  satisfy 

Uk  =  |j~jjj i/u(x)pk(x)w(x)dx.  (1.3.4) 

For  an  integer  n>0,  the  truncated  series  of  u  of  order  n  appears  as 

Pnu  =  SiikPk(x).  (1.3.5) 

k*0 

1.4  GaussLobatto  Quadratures  and  Discrete  Polynomial  Transforms 

Expanding  any  u(x)eL2(— 1,1)  in  terms  of  the  coefficients  Uk,  called  the 
continuous  expansion ,  depends  on  the  known  function  u(x).  With  u(x)  is  not 
known  a  priori,  a  discrete  expansion  for  u(x)— which  depends  on  the  values  at  the 
nodes-must  suffice. 

A  close  relation  exists  between  orthogonal  polynomials  and  Gauss-Lobatto 
integration  formulas  on  the  interval  [—1,1].  Let  Xo,-.  ,Xn  equal  the  roots  of  the 
(n+l)—th  orthogonal  polynomial  pn*i  and  let  w0,...,wn  equal  the  solution  of 
the  linear  system  given  by 

E  (xj)kwj  =  f  xkw(x)dx  0  <  k  <  n,  (1.4.1) 

j  i0  »  -! 

where  w(x)  equals  the  weight  function  associated  with  the  Sturm-Liouville 
problem,  Wj>0  for  j  =  0,...,n,  and 

j?0P(xj)wj  =  /_1P(x)w(x)dx,  (1.4.2) 

hold  for  all  peP2n»i-  The  positive  numbers  Wj  are  called  "weights”  (see  Canuto 
et.  al.(1988)  for  proof).  This  version  of  Gauss  integration  produces  roots, 


Identifying  the  function  u(x)  as  "square  integrable"  on  the  given  domain 
requires  /|u(x)J2dx  <  ». 
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corresponding  to  the  collocation  points,  which  appear  in  the  interior  of  (—1,1). 
Since  boundary  conditions  require  one  or  both  end  points,  a  generalized  Gauss 
integration  formula  must  include  these  points. 

The  Gauss-Lobatto  formula  considers 

q(x)  =  pn»i(x)  +  apn(x)  +  bpn-i(x),  (1.4.3) 

with  a  and  b  chosen  so  that  q(— l)=q(l)=0.  For  a  given  weight  function 
w(x)  and  corresponding  sequence  of  orthogonal  polynomials  pk,  k=0,  1,  2,... 
we  denote  by  x0,...,xn  the  nodes  of  the  n+1  point  integration  formula  of 
Gauss-Lobatto  type,  and  by  wo,...,wn  the  corresponding  weights. 

In  a  collocation  method  the  fundamental  representation  of  a  smooth 
function  u  on  (—1,1)  appear  in  terms  of  its  values  at  the  discrete  Gauss- 
Lobatto  points.  Approximate  derivatives  of  the  function  occur  by  analytic 
derivatives  of  the  interpolating  polynomial.  The  interpolating  polynomial, 
denoted  by  Inu,  belongs  in  the  set  Pn  and  satisfies 

Inu(xj)  =  u(xj)  0  <  j  <  n.  (1-4  4) 

Since  (1.4.4)  represents  a  polynomial  of  degree  n,  it  admits  an  expression  given 

by 

Inu  =  2  UkPk(x).  (1.4.5) 

k  *  0 

Since  the  interpolating  polynomial  must  satisfy  the  function  exactly  at  the  nodes, 
we  get 

u(xj)  =  k£oukpk(xj),  (1-4.6) 

where  Uk  equal  the  discrete  polynomial  or  expansion  coefficients  of  u.  The 
inverse  relationship  satisfies 

Uk  =  ^  j£Q  u(xj)pk(xj)wj,  (1.4.7) 
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where  the  coefficients  7k  equal 


7k=  p2  (xj)wj.  (1.4.8) 

Equations  (1.4.6)  and  (1.4.7)  enabling  transforms  between  physical  space 
{u(xj)}  and  spectral  space  {iik}  are  called  the  discrete  polynomial  transforms 
associated  with  the  weights  wo,...,wn  and  the  nodes  xo,...,xn. 

1.5  Legendre  Polynomials 

A  collection  of  the  essential  features  of  Legendre  polynomials  appears 
below.  The  Legendre  polynomials  {Lk(x),  k  =  0,  1,...,}  equal  the  eigenfunctions 
of  the  singular  Sturm-Liouviile  problem  given  by 

g^(l-x2)jj£Lk(x)]  +  k(k  +  l)Lk(x)  =  0,  (1.5.1) 

which  equals  (1.1.1)  with  p(x)=l— x2,  q(x)=0  and  w(x)=l.  By  normalizing 
Lk(x)  so  that  Lk(l)=l,  the  Legendre  polynomials  satisfy 

Lk(x)  =  jkjjri  ^rk  (x2  -  l)k  (1.5.2) 


and  represent  the  solution  to  (1.5.1)  with  boundary  conditions  (1.5.5).  These 
polynomials  also  satisfy  the  recurrence  relation,  expressed  as 


Lk»i(x)  —  k  j-  xLk(x)  jj.  j-  Lk-i(x), 

(1  5.3) 

where  Lo(x)=l 

and  Li(x)=x, 

|Lk(x)|<l,  -1  <  x<  1, 

(1.5,4) 

Lk(*l)  =  (±l)k, 

(1.5.5) 

|L'k(x)|  <  jk(k  +  1),  -1  <  x  <  1, 

(1.5.6) 

L/k(±l)  =  (*l)k  ^k(k  +  1), 

(1.5.7) 

and 

J*  |  L£(x)dx  =  (k  +  j)'1, 

(1.5.8) 

along  with  the  property  that  Lk(x)  is  even  if  k  is  even,  and  odd  if  k 

is  odd. 

The  continuous  expansion  of  any  ueL2(-l,l)  in  teims  of  the  Legendre 
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polynomials  appears  as 


u(x)  =  E  UkLk(x).  (1.5.9) 

Multiplying  both  sides  by  Lj(x)  and  integrating  from  x=-  1  to  x=l,  gives 

f  |  u(x)Lj(x)dx  =  j  i  UkLj(x)Lk(x)dx,  (1.5.10) 

where  w(x)=l  according  to  (1.3.2).  Using  the  orthogonal  properties  of  the 
Legendre  polynomials  and  (1.5.8)  gives 

uj  =  (j  +  j)  /^(x)Lj(x)dx.  (1.5.11) 

The  Legendre  polynomials  may  appear  directly  as  the  expansion  functions  of 
(1.5.9)  or  as  a  combination  of  Legendre  polynomials  which  satisfy 

V'j(xi)  =  %  (1-5.12) 


where  ip(x)  equals 

^  W  =  n(n-f  l)Ln(xj)  x1  -  *xj  3xL,,^x^  (1.5.13) 

This  later  form  allows  simpler  implementation.  The  nodes  {xj},  j=l,...,n— 1, 

equal  the  zeros  of  ^Ln(x),  with  x0=-  1,  and  xn=l.  The  quadrature  weights, 

shown  in  (1.4.2),  satisfy 

Jo^j(xk)wk  =  /j  ^(x)w(x)dx.  (1.5.14) 

Since  the  Legendre  polynomials  correspond  to  w(x)=l,  and  the  expansion 
polynomials  satisfy  the  relation  ^j(xk)=4j,  (1.5.14)  reduces  to 

k?0  ^kjwk  =  f  x  V'jWdx  (1.5.15) 

or 

wj  =  f  {  ^j(x)dx.  (1.5.16) 

Inserting  the  expression  for  the  expansion  polynomial,  (1.5.13),  gives  after 
integration 


-  ^2+  l)L.(i1)=  ■  J“  0 . »  <1!u7> 

for  the  quadrature  weights.  The  normalization  factors  7k,  introduced  in  (1.4.8) 
for  the  general  Sturm-Liouville  problem,  equal 

7k  =  (k  +  j)‘\  for  k  <  n  (1.5.18) 

and  7n  =  2/n.  (1.5.19) 

foT  the  specific  polynomials  given  in  (1.5.13). 

1.6  Differentiation  Using  Legendre  Polynomials 

Differentiation  may  occur  in  either  spectral  space  or  physical  space. 
Differentiation  i  u  .  •  .cral  space  consists  of  computing  the  Legendre  expansion  of 
the  derivative  of  a  function  in  terms  of  the  Legendre  expansion  of  the  function 
itself.  For  example,  if  u(x)=£k=o  UkLk(x),  ^  can  be  represented  as 

^  =  <5  =  k?o3^k)Lk(x)l  (L61) 

where 


j^fik  =  (2k  +  1)  up,  P+k  odd.  (1.6.2) 

The  proof  of  (1.6.2)  begins  with  a  relation  between  Legendre  polynomials  and 
their  derivatives: 

(2k  +  l)Lk(x)  =  ^Lk*i(x)  —  g^Lk-i(x),  k>0.  (1.6.3) 

Substituting  this  expression  for  Lk(x)  into  (1.6.1)  gives 


Jo  SF+HT  35EI'k*|(x) 

l  diik/dxd  T _  /_x 

_k?o5i+TaiLk*lw 

which  upon  changing  the  limits  of  the  summation  gives 

%  duk«  i  /dx  d  r  _f„\ 

_k?-i^+Va^Lkw- 

Combining  both  terms  gives 


(1.6.4) 


(1.6.5) 
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|.uW  =  k|[^dx_du^dx]d.Ll(l) 


(1.6.6) 


where  both  the  terms  corresponding  to  k=0  and  —1  vanishes  since  ^L0(x)=0 
and  ^L.j(x)=  0.|  This  expression  represents  the  derivative  of  u(x)  in  spectral 
space. 

In  physical  space,  the  derivative  appears  as 


Hx^Su)  =  dx1^  =  k?0fik3xLkM- 

Equating  (1.6.6)  and  (1.6.7),  and  recognizing  that  g^L0(x)=0,  gives 


CO 


CO 


k?i“kaacLkW=k51 


duk-i/dx  duk*  1  /dx 

51  +  3" 


3^Lk(x)- 


(1.6.7) 


(1.6.8) 


Since  the  Lk(x)  equal  the  eigenfunctions  of  Sturm-Liouville  problem,  they,  along 
with  their  derivatives,  form  a  linearly  independent  or  orthogonal  set.  Therefore, 
multiplying  (1.6.8)  by  ^Lj(x)  and  integrating  from  -1  to  1  results  in 


Uk  2k  -  1  2k  +  3  ’ 
Evaluating  (1.6.2)  with  n=k+l  gives 

(1.6.9) 

d  ® 

^un-i  *  (2n  -  1)  p?nup,  p+n  even; 

similarly,  (1.6.2)  evaluated  using  n=k— 1  yields 

(1.6.10) 

d  ® 

=  (2n  +  3)  E  up,  p+n  even. 

p=nf2 

Substituting  (1.6.10)  and  (1.6.11)  into  (1.6.9)  gives 

(1.6.11) 

0D  W 

iin  =  S  Up  -  E  Up,  p+n  even, 

p=n  p=n«2 

which  reduces  to  un  =  un.  This  completes  the  proof  of  (1.6.2). 

(1.6.12) 

The  two 

derivative  expressions  given  by  (1.6.1)  and  (1.6.7)  produce  different  results  in 
practice: 


$For  k=— 1  we  have  dL.j(x)/dx  =  c/(l — x2),  where  c  equals  a  constant.  Since 

the  boundary  conditions  are  dLk(±l)/dx=(±l)kk(k+l)/2  which  equals  zero  for 
k= — 1,  we  see  that  the  constant  c  must  equal  zero. 
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axC*-®)  *  P"-E-  (1-6.13) 

The  quantity  on  the  left  equals  the  Legendre  Galerkin  derivative.  The  error, 
g^{Pnu)— Pn-ig^  decays  spectrally  for  infinitely  smooth  solutions.  However,  for 
functions,  u,  with  finite  regularity  (not  infinitely  periodic)  this  difference  decays 
at  a  slower  rate  than  the  truncation  error  for  the  derivative  a^_Pn-4l-  Thus 
^(Pnu)  is  asymptotically  a  worse  approximation  to  ^  than  Pn-gj  (see 
Canuto  et.  al.  (1988)). 

1.7  Legendre  Derivatives  at  the  Nodes 

Approximate  differentiation  in  physical  space  occur  by  differentiating  the 
interpolation  Inu  (as  defined  in  (1.4.5))  and  evaluating  it  at  the  nodes.  This 
resulting  polynomial  of  degree  n—  1,  represented  as 

DnU  =  (L71) 

and  called  the  Legendre  collocation  derivative  of  u  relative  to  the  chosen  set  of 
nodes,  differs  from  the  Galerkin  derivative  ^(Pnu)  since  the  latter  depends  on 
the  continuous  coefficients  iik  and  the  former  on  the  discrete  coefficients  fik- 
One  method  for  obtaining  the  collocation  derivative,  involves  computing 
the  values  (Dnu)(xi),  (i  =  0,...,n)  from  the  values  u(xj),  (j  =  0,...n),  by 
employing  (1.4.7)  for  the  discrete  Legendre  coefficients  uj,  and  (1.6.2)  for  the 
discrete  derivative  coefficients  and  computing  (Dnu)i  from 

(Dnu)(xi)  =  J^(uk)t&k(xi).  (1.7.2) 

A  preferred  option  involves  the  collocation  derivative  at  the  nodes  through 
matrix  multiplication.  It  appears  as 

(Dnu)(xi)  =  ^Uk  I  x=x.  (1-7-3) 

for  i=0,...,n.  When  Dik=^^(*i),  (1.7.3)  equals 
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(1.7.4) 


(Dn«)(xi)  =  E  UkDik,  for  i  =  0 . n. 

k  =  0 


Using  (1.5.13)  for  ^\(x)  gives 

Ln(xi 


Du  = 


tn{xk)  Xi-Xk’  1  *  k' 

i  =  k  =  0.  (1.7.5) 

_(n_+_l)n  i  =  k  =  n. 

0,  otherwise. 

1.8  Integration  Using  Legendre  Polynomials 

Integration  in  transform  space  consists  of  computing  the  integral  of  the 

Legendre  expansion  of  a  function.  If  u(x)a  En  Uk^k(x),  the  integral  over  the 

k=0 

domain  x€[— 1,1]  equals 


n(x)dx  a  J ’  J  uk^(x)dx.  (1.8.1) 

Assuming  the  series  converges,  integration  and  summation  may  change  places, 
giving 

J*  ^  u(x)dx  a  iik  J  ^  ^(x)dx.  (1.8.2) 

Using  the  integral  of  the  expansion  function  according  to  (1.5.16)  gives 

J  i  u(x)dx  a  X  UkWk-  (1.8.3) 


CHAPTER  H 

SPECTRAL  ELEMENT  METHOD 

2.1  Introduction 

The  spectral-element  method,  a  variational  procedure  in  which  the 
approximating  functions  depend  on  representing  the  given  domain  as  a  collection 
of  simple  sub-domains,  differs  from  both  spectral  methods  and  finite-element 
methods  in  two  ways:  (1)  pure  spectral  methods  employ  high  degree 
approximating  functions  with  support  defined  over  the  entire  domain,  and  (2) 
finite-element  methods  use  low  degree  approximating  functions  with  compact 
support  (i.e.,  a  given  element’s  approximating  functions  differ  from  zero  only 
within  the  element).  Spectral-element  methods  exploit  the  advantage  of  high 
degree  functions  inherent  in  pure  spectral  methods,  along  with  the  flexibility 
finite-element  methods  provide  in  representing  complex  domains.  The  sub- 
domains,  or  finite  elements,  equal  geometrically  simple  shapes  that  permit  a 
systematic  construction  of  the  approximating  functions.  These  ecumenical 
functions  satisfy  all  boundary  conditions  and  problem  data  by  employing 
concepts  of  orthogonal  polynomials  from  Sturm-Liouville  theory.  On  an 
elemental  basis,  the  dependent  variables  appear  as  a  finite  sequence  of  the 
approximation  functions  with  coefficients  representing  the  dependent  variables  at 
a  finite  number  of  preselected  points  (i.e.,  nodes,  whose  number  and  location 
dictates  the  degree  and  form  of  the  approximating  functions). 

2.2  Partitioning  of  Domain 

One  feature  of  the  spectral-element  method  distinguishing  it  from  the  pure 
spectral  method  allows  representing  the  given  domain  by  a  collection  of  sub- 
domains.  A  subsequent  transformation  maps  each  sub-domain  from  the  physical 
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(x,y,z)  space  to  the  local  (r,s,t)  space  by  an  isoparametric  mapping.  The  sub- 
domains  in  local  space  equal  simple  geometries,  such  as  cubes  in  three- 
dimensional  space.  Two  important  features  in  typical  geometries  dictate  this 
mapping:  first,  the  definition  of  the  approximation  functions  from  Sturm- 
Liouville  equations  only  apply  to  certain  well-defined  geometries;  and  second,  an 
arbitrary  domain  cannot  accept  a  collection  of  simple  domains  without 
introducing  error.  By  defining  the  approximating  functions  element-wise,  the 
accuracy  of  the  approximation  improves  by  increasing  either  the  number  of 
elements  (i.e.,  refining  the  mesh)  or  the  degree  of  the  approximating  functions. 

In  mathematical  terms,  the  total  domain  n=flUcK)  splits  into  a  finite 
number,  E,  of  subsets,  Cle,  called  finite  elements,  such  that:  each  fie  is  closed 
and  non-empty;  the  boundary  dfte  of  each  Qe  is  Lipschitz-continuous  (no 
singularities,  cusps,  et  cetera);  the  intersection  of  any  two  distinct  elements  is 
empty,  i.e.,  ftenflf=®,  e#;  and  the  union  fi  of  all  dements  fie  equals  the 
total  domain,  given  as 

E 

n  =  Ene-  (2.2.1) 

e  =  1 

We  could  not  satisfy  the  last  property  without  the  mapping  between  physical  and 
local  space. 

2.3  Spectral-Element  Interpolation 

By  allowing  the  possibility  that  each  element  represents  the  entire  domain 
with  the  general  boundary  conditions  of  the  differential  equation,  the  essential 
boundary  conditions  equal  the  values  of  the  independent  variables  at  the  nodes, 
while  the  natural  boundary  conditions  get  subsumed  into  the  variational  form  of 
the  equation  over  the  element.  After  assembling  the  elements,  the  boundary 
values  on  portions  of  the  boundaries  of  elements  sharing  the  boundary  of  tho 
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given  domain  are  replaced  by  the  actual  specified  values  (imposition  of  boundary 
conditions). 

In  the  spectral-element  method,  the  minimum  degree  of  the  algebraic 
approximating  functions  depends  on  the  order  of  the  differential  equation  being 
solved,  and  the  degree  of  the  polynomial  in  turn  dictates  the  number  of 
interpolation  points,  called  nodes,  to  be  identified  in  the  element. 

The  approximation  functions,  also  called  interpolation  functions,  depend  on 
interpolation  of  the  function  and  possibly  its  derivatives  at  the  nodes  of  the 
element.  The  nodes  placed  along  the  boundary  of  the  element  uniquely  define  the 
element  geometry.  Place  any  additional  nodes  required  to  define  the 
interpolation  functions  at  other  points,  either  in  the  interior  or  on  the  boundary. 
The  boundary  nodes  also  enable  the  connection  of  adjacent  elements  by  requiring 
equality  of  the  primary  degrees  of  freedom  (i.e.,  variables  that  appear  in  essential 
boundary  condition)  at  nodes  shared  by  any  two  elements.  Thus,  we  cannot 
accurately  represent  discontinuous  primary  variables.  Such  problems  arise  in,  for 
example,  the  study  of  compressible  flow  where  shock  waves  contain  velocity 
discontinuities.  These  functions  make  poor  primary  variables  in  the  spectral- 
element  model  unless  we  employ  special  procedures  during  assembly. 

For  each  fte,  let  PS  denote  the  finite-dimensional  spaces  spanned  by 
linearly  independent  local  interpolation  functions  {^?}i=o  of  the  nodal  points. 
Over  each  element  fiecn  the  approximation  Inue  of  ue  equals 

u*  »  Inue  =  ,£u?^(r),  (2.3.1) 

where  the  local  expansion  coefficients  u?  equal  the  values  of  ue  at  the 
preselected  nodes  {r?}  in  the  element  Cle.  As  indicated  in  §1.5,  the 
interpolation  functions  satisfy 
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(2.3.2) 


$(fj)  =  «5ij, 

where  6jj  is  the  Kronecker-delta  function. 

2.4  Connectivity  (or  Assembly)  of  Elements 

As  mentioned  earlier,  all  elements  contain  boundary  nodes  defining  their 
geometry  and  allowing  connection  with  their  neighbors.  The  connectivity  of 
elements  requires  equal  values  of  the  primary  variables  in  nodes  common  to 
adjacent  elements.  Assembling  the  unique  sub-domains  into  the  entire  domain,  a 
process  known  as  direct  stiffness ,  requires  the  identification  of  a  universal  or 
global  system  of  nodes,  and  a  corresponding  set  of  global  expansion  coefficients 
for  the  primary  variables.  The  resulting  matrix  expression  relates  the  global 
expansion  coefficients  to  the  parameters  of  the  governing  differential  equation 
and  boundary  conditions. 

As  indicated  earlier,  the  primary  variables,  those  associated  with  the 
essential  boundary  conditions,  appear  as  the  global  expansion  coefficients  of  the 
assembled  matrix  relation.  The  secondary  variables  appear  in  the  natural 
boundary  conditions. 

2.5  Isoparametric  Formulation 


Isoparametric  schemes  use  the  same  interpolation  functions  to  represent 
both  the  co-ordinate  mapping  and  the  primary  variables.  Thus,  the  physical 
space  x  maps  into  the  local  r  co-ordinate  system  by 


xe  s  Inxe  =  J0*W(r), 

(2.5.1) 

when  primary  variables  appear  as 

ue  s  Inue  =  E  u?^?(r). 

is  0 

(2.5.2) 
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CHAPTER  m 

TIME-SPLITTING  SCHEME  FOR  THE  NAVIER-STOKES  EQUATIONS 
3.1  Governing  Equations  and  Boundary  Conditions 

According  to  the  formulation  given  below,  the  domain  under  consideration, 
flcR3  with  boundary  dfi,  may  translate  but  not  deform.  The  Navier-Stokes 
equations  on  the  closed  domain,  n=QUdft,  consist  of  the  constant  density 


continuity  equation  given  by 

v*v  =  o, 

(3.1.1) 

and  the  corresponding  momentum  equation  expressed  as 

^+(V*V)V  =  -^+f  +  i/V*V. 

(3.1.2) 

Introducing 

V*V*V  =  -<V-V)V  +  ^(V-V) 

(3.1.3) 

into  the  momentum  equation  gives 

j~  =  V*V*V  -  V(£  +  iV*  V)  +  f  +  I/V2V, 

(3.1.4) 

in  which  both  V  and  p  along  with  the  body  force,  f,  depend  on  both  position 
in  the  fluid  and  time. 


The  physical  boundary  conditions  for  a  given  problem  must  allow  us  to 
divide  the  entire  boundary  into  regions  associated  with  essential  boundary 
conditions  (e.g.,  walls  and  inlets),  natural  boundary  conditions  (e.g.,  outlets  and 
free-streams),  and  per  odicity  boundary  conditions.  The  numerical  procedure 
handles  each  of  these  regions  separately. 

The  essential  boundary  conditions  appear  as 

V  =  VwaU  (3.1.5) 

on  solid  wall  boundaries,  and 

V=Vin  (3.1.6) 
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on  inlet  boundaries.  The  natural  boundary  condition  at  the  outflow  equals 

(n-V)V  =  0,  (3.1.7) 

where  n  is  the  outward  surface  normal  vector,  and  along  the  free  stream 

n*V  —  n-Vboundary  (3.1.8) 

and 

n*(£-n)  =  0,  (3.1.9) 

where  a  equals  the  local  stress  tensor.  Equation  (3.1.9)  expresses  the  fact  that 
the  stress  at  the  free  stream  must  lie  entirely  normal  to  the  boundary.  By 
manipulating  the  pair  of  conditions  (3.1.8)  and  (3.1.9)  we  can  show  that  they 
equal  the  homogeneous  natural  condition  on  velocity,  (n*V)V=0,  (see  appendix 
C).  The  boundary  condition  along  periodic  surfaces  equals 

V(x+L)  =  V(x),  (3.1.10) 

where  L  equals  the  relative  position  vector  between  the  two  periodic  boundaries 
All  of  these  conditions  refer  to  velocity;  pressure  boundary  conditions  depend  on 
the  governing  equations.  Finally,  the  initial  conditions  equal  V(x,t=0)=V0(x) 
for  xeflcR3. 

3.2  Splitting  Method 

In  the  variational  solution  of  time-dependent  problems,  we  represent  the 
dependent  variables  in  a  finite  dimensional  vector  space.  The  undetermined 
coefficients  depend  on  time,  while  the  base  functions  depend  on  spatial  co¬ 
ordinates.  This  leads  to  a  two-stage  approximation,  both  of  which  could  employ 
variational  methods.  We  choose  to  discretize  the  equations  in  space  and  iterate 
in  time,  thus  giving  a  spatial  variational  problem.  Such  a  procedure, commonly 
known  as  a  semi-discrete  approximation,  results  in  a  set  of  ordinary  differential 
equations  in  time. 
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No  variational  form  exists  for  the  full  momentum  equation;  therefore,  we 


split  it  into  simple  forms  and  apply  variational  techniques  to  each  portion 
individually.  Employing  the  splitting  scheme  followed  by  Wessel  (1992),  we 
introduce  intermediate  velocities,  V,  V,  V,  and  V,  which  allows  splitting  of 
the  momentum  equation  into  fractional  steps.  The  scheme  employs  a  "predictor- 
corrector"  approach  whereby  the  predicted  velocity  at  time  step  n+1,  which 
results  from  sequentially  computing  intermediate  velocities  based  on  the  non¬ 
linear  terms,  the  viscous  terms,  and  the  pressure  terms,  determines  the  pressure 
gradient.  The  corrected  velocity  depends  on  this  predicted  pressure  gradient.  A 
brief  explanation  follows.  The  first,  or  non-linear,  step  appears  as 

Vn  +  l  _  \ra 

■■■  —  =  VM*V»  +  F+1;  (3.2.1) 

the  second,  or  viscous,  step  equals 


yiMl  _  yn  *1 
At 


=  1/V2V”; 


(3.2.2) 


and  the  third,  or  pressure,  step  includes 

—  -t  -  =  -V(£  +  =  -VPnM.  (3.2.3) 

The  velocity  after  the  pressure  step,  V71*1,  must  satisfy  the  divergence-free 
constraint.  By  applying  the  divergence  operator  to  (3.2.3),  a  relationship 
between  pressure  and  V”*1  ensues.  The  solution  of  this  Poisson  equation  for 
pressure  determines  the  predicted  pressure  gradient.  Using  the  previously 
computed  velocity  after  the  non-linear  step  and  this  new  pressure  gradient,  a  new 
velocity  results  after  adding  the  explicit  viscous  term: 


— - —  =  -VPn  ♦>  +  M2V".  (3.2.4) 

We  still  must  add  another  to  the  right-hand  side  to  yield  the  full  Navier- 

Stokes  equation.  We  use  the  Crank-Nicholson  method  by  adding  to 
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the  velocity  V0  *l  giving  an  implicit  viscous  step  appearing  as 

=  jrfiV**1.  (3.2.5) 

We  benefit  from  this  splitting  scheme  by  representing  the  implicit  viscous  and 
pressure  steps  in  time-independent ,  elliptic  form:  the  viscous  step  in  the  form  of 
Helmholtz’s  equation,  and  the  pressure  step  as  Poisson’s  equation.  We  express 
both  equations  in  variational  form  and  solve  them  by  finding  the  extremum  of 
the  corresponding  functional.  Some  of  the  details  of  each  of  the  five  steps  appear 
below.  The  two  steps  containing  the  explicit  viscous  terms  require  no  exposition. 

3.3  Non-Linear  Step 

We  solve  the  non-linear  advective  term  explicitly  using  a  three-step  Adams- 

Bashforth  method  given  by  (3.3.1) 

=  b0(v*v*v  +  f)n+  Bj(v*v*v  +  f)n-1+  b2(v*v*v  +  f)n~24 

This  hyperbolic  operator  imposes  stability  conditions,  in  the  form  of  a  Courant- 
Friedrich-Lewy  number,  on  the  scheme.  Neither  bour  ^arv  conditions  nor 
continuity  constraints  apply  to  V11*1. 

3.4  Pressure  Step 

The  velocity  after  the  pressure  :tep,  V,  must  satisfy  the  zero  divergence 
constraint.  Applying  the  divergence  operator  to  (3.2.2)  gives 

=  V-(-VPnn).  (3.4.1) 

Since  V  does  not  satisfy  the  zero  divergence  constraint,  (3.4.1)  simplifies  to 

=  V*VPntl.  (3.4.2) 

u  V 

The  boundary  conditions  on  pressure  depend  on  the  governing  equations. 


$The  three-step  Adams-Bashforth  coefficients  equal  P0=23/12,  B j=— 16/12, 
and  B2=5/12. 
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We  obtain  the  pressure  boundary  conditions  by  taking  the  inner  product  of 
(3.1.4)  with  the  surface  outward  unit  normal  n  and  rearranging.  This  yields 
n-VP  =  — ^{n- V)  +  im-W  +  n-f  +  n-(V*V»V).  (3.4.3) 

This  expression  represents  the  exact  physical  boundary  condition  on  pressure 
obtained  from  the  governing  differential  equation.  Unfortunately,  we  cannot  use 
this  form  since  some  of  the  terms  on  the  right-hand  side  remain  unknown. 
Therefore,  the  numerical  procedure  uses  a  simplified  form  which  neglects  the 
viscous  term,  the  body  force,  and  the  non-linear  term: 

n-VP  =  -|-(n-V).  (3.4.4) 

When  we  discretize  this  equation  in  time  by  writing  the  time  derivative  of 
velocity  as  (V1**1— V^/nt,  we  see  that  unknown  terms  still  remain  on  the  right- 
hand  side.  By  further  approximating  the  time-derivative  using  intermediate 
velocities  we  obtain 

n-VP  =  -n.(Vnn-VD*1)/At.  (3.4.5) 

It  appears  that  we  still  do  not  know  V11*1  at  this  point;  however,  a  careful 
vetting  of  the  various  boundary  conditions  reveals  otherwise.  This  mathematical 
approximation  to  the  physical  boundary  condition  may  be  good  or  bad  depending 
on  the  characteristics  of  the  flow  as  indicated  below.  Where  essential  boundary 
conditions  occur  n-V11*1  equals  n*Vwaii  or  n-Viniet,  and  (3.4.5)  becomes 

n-VP  =  — n*(Vwaii  -  V  (3.4.6) 

This  term  differs  from  the  true  physical  boundary  condition  at  a  solid  surface  (or 
at  the  inlet  if  n-(V*V*V)=0)  since  it  lacks  both  the  viscous  term,  im-y2V,  and 
the  body-force  term,  n*f.  For  a  flow  with  no  body  force,  the  approximate 
boundary  condition  differs  from  the  physical  boundary  condition  by  the  viscous 
term.  In  other  words,  the  mathematical  boundary  condition  equals  the  inviscid 
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boundary  condition  along  solid  walls.  Deville  and  Orszag  (1980)  showed  that 
this  approximation  introduces  a  time-splitting  error  0(1)  in  n-V2V  over  a 
layer  of  thickness  0(<J& tv).  No  error  estimate  exists  at  the  outlet  since  even  the 
physical  boundary  conditions  remain  unknown. 

Along  the  free-stream  the  viscous  term  remains  negligibly  small  under  most 
circumstances;!  while  the  advection  term,  n*(V*¥*V),  equals  zero,  since  n 
and  V*V*V  are  perpendicular. %\  Furthermore,  when  the  body  force ’s 
perpendicular  to  the  free-stream  boundary— as  in  most  flows— the  body-force  term 
vanishes,  and  the  simplified  boundary  condition  equals  the  physical  one. 

Along  exit  boundaries  where  natural  conditions  get  specified,  the  terms 
neglected  by  the  approximate  pressure  boundary  condition  do  not  vanish. 
However,  their  influence  remains  confined  to  a  region  near  the  boundary,  since 
convection  sweeps  any  induced  errors  out  of  the  domain. 

3.5  Viscous  Step 

The  implicit  viscous  step  appears  as 


|The  viscous  term  equals  n»V2V.  The  velocity  can  be  written  in  terms  of  a  local 
co-ordinate  system  where  n  equals  the  normal  to  the  surface  and  s  lies 
tangential  to  the  surface.  Hence  V=Vss-fVnn,  and  n-V2V= 
n«[(d2Vs/5s2+d2Vs/dn2)s+(52Vn/9s2+d2Vn/dnJ)n].  The  first  term  is  zero 
because  s  and  n  remain  perpendicular.  The  second  term  vanishes  in  most 
common  flows.  For  example,  when  a  homogeneous  body  force  or  a  body  force 
with  a  vanishing  or  zero  gradient  near  the  free-stream  boundary  (e.g.,  the 
"Blasius"  body  force)  forces  the  flow,  the  second  derivatives  of  Vn  vanish 
sufficiently  far  from  disturbance  generating  structures.  Likewise,  when  the  free- 
stream  boundary  represents  a  moving  plate,  as  in  Couette  flow,  d2Vn/ds2 
usually  equals  zero,  and  d2Vn/dn2  again  vanishes  far  away  from  the  disturbance 
generating  structures. 

JfThe  boundary  conditions  specified  along  a  free  stream  reduce  to  n*VV=0, 
according  to  the  results  of  appendix  C.  Thus,  the  gradient  of  velocity  lies 
parallel  to  the  surface,  or  alternatively,  the  cross  product  of  V  and  V  lies 
parallel  to  n  (the  zero  normal  velocity  constraint  requires  V  to  lie  in  the  plane 
of  the  boundary).  Therefore,  the  cross  product  of  V  and  V*V  lies  in  the  plane 
of  the  surface,  or  perpendicular  to  n.  Whence  n>(V*V*V)=0. 
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=  *iA?2V®41. 


V® 


♦  1 


yn  < 


At 


(3.5.1) 


This  implicit  equation  remains  unconditionally  stable.  Therefore,  we  could  avoid 
unreasonable  time  step  restrictions  due  to  the  high  spatial  resolution  of  spectral 
approximations  near  the  boundaries  of  elements  save  the  other  explicit  steps. 

The  boundary  conditions  on  the  velocity  after  the  viscous  step,  V®41,  equal  the 
physical  boundary  conditions  given  in  §3.1.  In  formulating  the  viscous  step,  we 
did  not  subsume  the  zero  divergence  condition  into  the  expression  for  V®*1. 
Hence,  V®*1  does  not  satisfy  the  Navier-Stokes  equations  exactly.  Normally, 
this  error  does  not  dominate  the  solution  since  the  velocity  after  the  viscous  step 
nearly  equals  the  zero  divergence  velocity,  V11*1.  Amon  (1988)  observed  that 
the  divergence  of  V™1  remains  a  few  orders-of-magnitude  smaller  than  that  of 
V®41.  We  transform  (3.5.1)  into  Helmholtz’s  equation  by  adding  —  V®*1  to 
both  sides  of  (3.5.1)  giving 

-  V®41  =  -  V®41  +  *i/AtV*V®41.  (3.5.3) 

Rearranging  yields 

<3-5-4) 

which  appears  as 

-  Aty  =  F,  (3.5.5) 

when  <|>=V®41,  A2=^  and  F=  -A2V®41.  Since  $  depends  on  the  velocity 
satisfying  the  physical  boundary  conditions,  the  condition  on  <j>  associated  with 
essential  boundary  conditions  equals 

$  =  Vwall  or  V inlet-  (3.5.6) 

The  exit,  or  natural,  boundary  condition  equals 

n*V$  =  0,  (3.5.7) 

which  also  satisfies  the  free-stream  boundary  conditions  when  the  free-stream 
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condition  on  velocity  behaves  according  to  the  restrictions  given  in  appendix  C. 
The  periodic  boundary  conditions  equal  <J>(x+L)=4(x),  where  L  equals  the 
relative  position  vector  between  the  two  boundaries. 
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CHAPTER  IV 
NON-LINEAR  STEP 

4.1  Introduction 

The  non-linear  step,  the  only  explicit  part  of  the  three-step  time-splitting 
scheme,  introduces  a  time-step  size  stability  restriction.  Since  only  first-order 
derivatives  appear,  no  benefit  accrues  from  casting  the  governing  equation  in 
"weak"  form;  instead,  we  apply  a  collocation  variational  approach.  (When  the 
governing  equation  has  an  equivalent  “weak"  form,  we  may  multiply  it  by  a 
suitably  differentiable  test  function,  integrate  over  the  appropriate  domain,  and 
then  integrate  by  parts.  The  resulting  variational  form  contains  derivatives  of 
lower  order  than  the  original  equation.)  In  the  collocation  technique,  the  test 
function  equals  the  Dirac-delta  function.  Since  this  function  has  no  derivative, 
the  resulting  solution  must  contain  as  many  derivatives  as  the  order  of  the 
governing  differential  equation.  For  this  first-order  equation,  therefore,  the 
solution  occupies  the  space  of  functions  H^fl). 

4.2  Variational  Form 

The  Adams-Bashforth  three-step  procedure  applied  to  the  non-linear 
portion  of  the  split  Navier-Stokes  equation  appears  in  §3.3  and  equals 

V**1  -V”  =  BoM(V*V*V  +  f)n+  B,at(V«V*V  +  f)n-‘+  B2At(V*V*V  +  f)n'2. 

(4.2.1) 

To  transform  this  equation  into  variational  form,  we  multiply  by  a  test  function, 
V^H(fl),  and  integrating  over  the  entire  domain,  giving 

f  [V11*1  -V11  -  B&At(V*V*V  +  f)n-  B1at(V*7«V  +  f)11'1-  B2At(V*V*V  + 

f)"'2]^  =  0.  (4.2.2) 

The  collocation  method  uses  xi),  where  the  Xi’s  equal  the  locations  of 
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the  nodal  points  where  the  differential  equation  is  satisfied  exactly.  Introducing 
i>i=S{x—Xi)  into  (4.2.2)  and  integrating  results  in  the  discrete  form  of  the 


differential  equation  in  terms  of  the  unknown  expansion  coefficients  V^k  and 
Vijk-  Products,  such  as  V*V*V,  also  appear  in  terms  of  their  coefficients  at  the 
nodes.  The  result  equals 


V?Ji^V?jk=  (4.2.3) 

B0&t(V*\f»V  +  f)^k  +  B1at(V«V«V  +  ftfji  +  B*t(V«V«V  + 

We  explore  the  spatial  discritization  of  the  various  quantities  in  the  following 
sections. 

4.3  Spatial  Discritization  of  Vortirity 

Express  the  vortirity,  V*V,  in  Cartesian  co-ordinates  as 

»»v=  (If  -  +  (It  ■ -  If  +  (g  • -  .  (4-3- 0 

where  the  partial  derivatives  refer  to  the  global  co-ordinate  system.  We  must 
transform  these  derivatives  from  the  global,  or  physical,  co-ordinate  system  to 
the  local  (r,s,t)  system.  According  to  appendix  A,  the  partial  derivatives 
expressed  in  terms  of  the  local  co-ordinates  equal 


d  e 


_  red  e  ,  e 


+  s|5r  +te 


r*3F  +Sr3s 


d  e 


—  re£_e  J_  efd_  +e 


+sy^s  +ty^ 


and 


d  e  _e3  e  ,  c 

=  fVTC  +  8 
a  e  a  e 

=  r%  +s! 


*3t  ’ 
d  e 


(The  superscript  e  refers  to  the  element  number.)  Introducing 
into  (4.3.1)  gives 


(4.3.2a) 

(4.3.2b) 

(4.3.2c) 
these  expressions 

(4.3.3) 


JThe  time  step  number  replaces  the  element  superscript  for  convenience  in 
notation;  remember,  however,  that  whenever  a  subscript,  such  as  ijk,  represents 
the  node  number  in  three  dimensions  a  corresponding  superscript  should  indicate 
the  element  number  under  consideration. 
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yxV*_  rre*We  .  ce^6  ,  te^We  re5VC  se^Ve 
VxV  -  lryjF  +  S>TT  +  t>TT  ~  rW  Sjtf? 


«sv 


,,r«au'.  ,e«ue,  ,,9U'  -«W'  e«w«  eaw'v 

+1**37  +S^s_  +  tzr3T  “rxFT  Sx?T  ~txFr^ey 

,,re*Ve  ,  se£Ve  .  te*ve  re^Ue  __  Re^Ue  _  ,edVe, 

+(r*7F  +  Sx7T  +  tx7T  r>7F  ®y??  t3^T  ^6z- 

The  three  components  of  the  vorticity  equal 

/*e_  e^We  ,  -iWe  te5We_  Te0Ve_  se0Ve_  ,etfVe  ,4  ,  4a^ 

**-  ryjr  +  sy^T  +  l>rr  r*37  Sztt  tz?r »  (4.3.4a) 

<*€_  e5Ue ,  e^Ue ,  teaue  Tc^We_  e^We  .e5We  M 

<r-  IxdT  +  s W  +  Wyr  r*3T  SrFT  txFT  ’  (4.3.4b) 

and  fl-  re?Ve.  Ee^Ve  .  te^Ve  -egU6  3effJe_,  te^e  f4  3  4c) 

ana  U-  rx^7  +  Sx^j  +  lyrgj-  sy^j  '■yj^  ■  (4.J.4CJ 

The  x-vorticity  evaluated  at  the  local  co-ordinate  rj,Sj,tk  appears  as 

(Cx)ijk  =  (ry)ijk^-J-ijk  +  (sy)ijk|^rijk  +  (ty)ijk^“ijk 

—  (rz)ijkjjjr^7ijk  —  (Sz)ijk|rjjrijk  “  (tz)ijk^J-ijk-  (4.3.5) 

The  partial  derivative  with  respect  to  the  local  co-ordinate,  r,  equals  (see 
appendix  A) 

n^ijk  =  Dia^jb^kc-  (4.3.6) 

Similar  expressions  exist  for  the  s  and  t  components.  Substituting  into  (4.3.5) 

6*ves  (4.3.7) 

(Cx)ijk  =  £  [(ry)fjkDia<5jb^kc  4"  (sy)ijk^iaDjb^kc  4-  (tyJijk^a^jbDkcjVVijk 
abc 

?  [(rz)ijkDia^jb^kc  4"  (Sz)ijk^iaDjb^kc  4"  (tz)fjk^ia^jbf^kc]^ijk- 
abc 

Expressing  the  other  vorticity  components  this  way  gives  an  equation  for  the 
vorticity  vector  at  the  local  node  (ri,Sj,tk):  ^4  3 

£1jk  =  abc|  ^[(ry)ijkDia^jb^kc  4"  (sy)ijk^iaDjb^kc  4-  (ty)fjk^ia^jbl^kc]Wijk 
[(rz)ijkDia^jb^kc  4-  (Sz)fjk^ial^jb^c  4-  (tz)ijk^ia^jbDkc]^ijk|ex 
4"  ^[(rz)ijkDia^jb^kc  4"  (Sz)ijk^iaDjb^kc  4"  (tz)ijk^ia^jbDkc]Uijk 
~  [(rx)ijkDia^jb^kc  4-  (Sx)ijk^iaDjb^kc  4-  (tx)fjk^ia^jbDkc]Wijk|  ey 
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+  ^[(rx)  ljkDia^jb^kc  +  (sx)ijk^iaDjb^kc  +  (tx)ijk£ia^jbDkc]Vijk 
~  [(ry)ijkDia^jb^kc  +  (Sy)ijk^iaDjb^kc  +  (ty)fjk^ia£jbBkc]Ufjkj  Ozj. 
4.4  Spatial  Discritusation  of  Cross-Product  of  Velocity  with  Vortidty 
The  cross-product  of  the  velocity  with  the  vortidty  appears  as 

V«<=  (VC*-WCy)ex  +  (WCr-UCa)ey  +  (UCjr-V<*)e*.  (4.4.1) 

Using  the  discrete  expressions  for  the  velodty  and  vortidty  gives 

(V*()?jk  =  jjc{  (4.4.2) 

f  [[(fi)  ij  kDia^j  b^kc  +  (Sx)ijk£iaDjb^c  +  (tx)ijktfiafjb^kc]Vijk 
~  I(ry)ijkDia^jb^kc  +  (sy)ijk^iaUjb^kc  +  (ty)ijk^ia^jbBkc]Uijk  V*jk 
—  ^[(rz)  ijkUia^jb^kc  +  (Sz)ljk^iaUjb^kc  +  (tz)ijk^ia^jbDkc]Uijk 
—  [(rx) ljkDia^jb^kc  +  (Sj)ijk^iaUjb^kc  (tx)?jk^ia^jbUkc]Wfjkl  Wijkl 


Ox 


+  [[[(ry)?jkDia5jb*kC  +  (Sy)fjk^iaUjb^kc  +  (ty)ijk^ia^jbUkc]Wfjk 
“  I(rz)ijkDia^jb^kc  +  (Sz)tjk^iaUjb^kc  +  (tz)ijk^ia^jbDkc]Vijk|  Wijk 

~  ^[(rx)?jkUia<5jb(5kc  +  (sx)ijk^iaUjb^kc  +  (tx)ijk^ia^jbUkc]Vfjk 
~  [(ry)ijkUia5jb^kc  +  (Sy)ljk^iaUjb^kc  +  ( ty )  i j k ^ia^j bUkc] Uf j k  Uijk  6y 
4"  J^|[(rz)ijkDia^jbfttc  +  (Sz)ijk^iaUjb^kc  +  (tz)ijk^ia^jbUkc]Uijk 

—  [(rx)1jkDia^jb^kc  +  (sx)ijk^iaUjb^kc  +  (tx)ijk^ia^jbDkc]Wfjkj  Ufjk 

—  [((ry)ijkD  ia^jb^kc  +  (s  y)  i j  k  <^iaD j  b  <^kc  +  (ty)fjk^ia^jb®kc]^ijk 
~  [(rz)ijkUia5jb^tc  +  (sz)ijk^iaDjb<5kc  +  (tz)ijk^ia^jbUkc]Vtjkj  Vjjkj  Czj 

for  the  cross-product  of  the  velocity  with  the  vortidty  at  the  local  node  (ri,Sj,tk). 

4.5  Summary 

The  equation  governing  the  velodty  after  the  pressure  step  equals 

=  V?jk  +  BoM[(V«C)?jk  +  ffjk]n  (4.5.1) 

+  B,At[(V*C)?jk  +  f!jk]n'3  +  B2at[(V«C)!jk  +  f!jk]n'2; 
where  the  terms  (V*£)ijk  appear  in  (4.4.2).  This  explicit  relation  for  the 
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unknown  velocity  coefficients  at  time  step  n+1  depends  on  known  quantities 
from  the  previous  step,  denoted  by  the  superscript  n.  Its  solution  does  not 
involve  inverting  the  stiffness  matrix;  unfortunately,  however,  this  simplicity 
comes  with  a  concomitant  loss  of  accuracy— the  collocation  scheme  searches  for 
the  solution  with  a  measurement  device  tolerating  first-order  errors  in  time. 
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CHAPTER  V 


PRESSURE  STEP 

5.1  Variational  Approach 

According  to  §3.4,  the  governing  differential  relation  for  the  pressure  step 
equals 

£|£-MP.  (5.1.1) 

This  represents  an  elliptic,  Poisson  equation;  consequently,  a  corresponding 
variational  form  exists.  For  the  moment,  assume  that  the  functional  whose 
Euler-Lagrange  equation  yields  (5.1.1)  equals 

J[P]  =  [-J  (VP-VP)at  +  V^-VPJdV.  (5.1.2) 

The  standard  boundary  conditions  accompanying  a  functional  of  this  type  equal 

n'WF=0,  (5.1.3) 

where  F  equals  the  integrand  of  the  functional.  Applying  this  boundary 
condition  formula  to  the  functional  shown  in  (5.1.2)  yields 

—  n-VPat  +  n-V11*1  =  0.  (5.1.4) 

According  to  the  derivation  in  §3.4  the  appropriate  boundary  conditions  for  the 
pressure  step  equal 

-  n-VPat  +  n-V71*1  =  n-V71*1,  (5.1.5) 

which  does  not  correspond  to  that  shown  in  (5.1.4),  since  an  additional  n-V11*1 
exists  on  the  right-hand  side;  therefore,  we  must  add  a  boundary  integral  to  our 
functional.  Thus,  the  functional  satisfying  the  appropriate  boundary  conditions 
and  the  governing  pressure  step  relation  equals  (5  16) 

J[P]  =  [-i  (VP-VP)at  +  V^-VPJdV  -  fdU  Ptn-V^da. 

Applying  the  Euler-Lagrange  equation  to  this  functional  yields  (5.1.1)  with 
boundary  conditions  (5.1.5).  A  demonstration  of  this  follows. 
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Up  until  this  point,  we  accepted  the  validity  of  the  functional  shown  in 
(5.1.6)  prima  facie.  We  must  show  that  the  extremum  of  this  functional  does 
indeed  yield  the  differential  relation  and  boundary  conditions  for  the  pressure 
step.  We  begin  by  considering  the  variation  of  J[P]  with  respect  to  P: 

fl[P]  =  Jfj  [-VP-V($P)At  +  Vn*1-V(^P)]dV  -  ^  (a-V^JAPdcr. 
Rearranging  gives  (5.1.7) 

<5J[P]  =  J rQ  [-V.(VP5P)at  +  flPV-VPAt  +  V.(V-WiP)  -  «PV-VI1+I]dV 

(5.1.8) 

Employing  Gauss’  divergence  theorem,  we  transform  the  first  and  third  volume 
integrals  into  surface  integrals  giving 

&J[P]  =  J^[V*VPAt  — V-V^^PdV  (5.1.9) 

+  fdQ  [-n-VPAt  +  n* V®*1  —  n*  V^tfPda. 

An  extremum  of  J[P]  occurs  when  this  variation  equals  zero.  Since  the 
variation  of  P  on  the  boundary  of  the  domain  does  not  depend  on  its  variation 
within  the  boundary,  the  extremum  occurs  when  both 

V-VPAt-V-V**1^,  (5.1.10) 

within  the  volume  of  the  domain,  and 

— n*VPAt  +  n-V"*1  -  n-V11*1  =  0,  (5.1.11) 

on  the  boundary,  are  satisfied  simultaneously.  We  see  that  these  equations  equal 
the  governing  equation  and  boundary  conditions.  Thus,  the  variation  of  the 
functional  given  by  (5.1.6)  yields  the  Poisson  equation  shown  in  (5.1.1)  with 
boundary  conditions  (5.1.5). 

5.2  Representation  in  Local  Co-ordinates 

We  must  transform  the  functional  in  (5.1.6)  from  the  global  x,y,z  co¬ 
ordinate  system  to  the  local  r,s,t  system.  The  corresponding  functional 
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representation  for  a  single  element,  designated  by  the  superscript  e,  in  local  co¬ 
ordinates  equals 

J[P]e  =  f  lJ  lJ  jH  (VP  *  VP)At  +  V-  VP]e  |  det[Je]  j  drdsdt 

-if'/-  j  P^n-^dAP,  (5.2.1) 

where  we  drop  the  superscript  indicating  the  time  step  number  and  replace  it  by 

a  superscript  indicating  element  number.  (Consult  appendix  B  for  details  of  the 

volume  element  conversion  from  global  to  local  co-ordinates.)  The  surface 

integral  contains  expressions  for  the  summation  over  the  six  faces  of  each 

element.  We  will  consider  the  details  of  these  terms  later. 

This  functional,  approximated  by  expressing  it  in  terms  of  the  expansion 

polynomials,  for  example,  Pea£  Ptbc^a(r)^b(s)^c(t),  appears  as 

a  be 

J[ple « f/'/l  i  V abcP abc '  V IbcP abcA t  (5.2.2) 

+  Vibe- VIbCPlbc]V’a(r)Vb(8)^c(t)  |  det[J«]  | abctirdsdt  -1  Pe. 

A  brief  explanation  of  the  accompanying  notation  seems  appropriate.  The 
subscript  abc  refers  to  the  node  ra,Sb,tc,  while  the  tilde  refers  to  the 

corresponding  unknown  expansion  coefficients.  We  postpone  treatment  of  the 

6 

surface  integral,  conveniently  expressed  in  the  interim  as  E  Pe,  until  §5.3. 

Many  of  the  following  numerical  minutiae  receive  a  more  thorough 
treatment  in  the  following  chapter  dealing  with  the  viscous  step.  The  gradient 
operator  in  global  co-ordinates  equals 

VIbcPIbc  =  (5  2.3) 

In  the  local  co-ordinate  system  this  operator  equals 

VIbcPIbc  =  (5-2.4) 

[dfl*bc(rx)abc  +  ^gbC(Sx)abc  +  "^^(^abclGx  + 
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[^rDabc  +  ^(s5)ab=  +  ^(tJUtcle,  + 

[^r^r*)abc  "*■  ^gbp(6»)abc  +  -^^{ti)abc]e*- 

Introducing  expansions  for  the  local  partial  derivatives  according  to  (6.3.5) 
through  (6.3.7)  gives 

VfbcPIbc  =  (5.2.5) 

,E^|[(rx)abc^ai^bj^ck  +  (Sx)abc^ail^bj^ck  +  (t*)abc^ai^bjI^ci]Cx  + 
((ry)abcl^ai^bj^ck  +  (Sy)abc^aiDbj  ^ck  +  (ty)abc^ai^bj®ck]®y  + 
[(rz)abcDai^bj  ^ck  +  (s|)abc^aiDbj^ck  +  (t|)abc^ai^bjDek]6zj  P?jk- 
The  inner  product  of  this  term  with  a  term  of  similar  form  appears  as 


^abcPabc*  VfbcPabc  =  (5.2.6) 

|[(rx)abcDai<5bj  <^ck  +  (sl)abc^aiDbj <^ck  +  (tx)abc<^ai<5bjDck]ex  + 
[(ry)abcDai^bj^ck  +  (sy)abc^aiDbj  ^ck  +  (ty)abc^ai^bjDck]6y  + 
[(rz)abcDai  5bj5ck  +  (si)abc^aiDbj  ^ck  +  (tz)abc^ai^bjDck]6z|  Pijk‘ 
lmn  [t(r*^a^>c®al^)ln  <5cn  +  (sx)abc^alDbm^cn  +  (tx)abc^al^b«Dcn]Cx  + 
[(ry)abcDal^bni<5cn  +  (Sy)abc^alDba^cn  +  (ty)abc^al^bnDcn]ey  + 
[(rl)abcDal^bm^cn  +  (Sz)abc^alDbB^cn  +  (tz)abc^al^bmDcn  Pfmn- 
Performing  the  inner  product  gives 


7|bcPIbc'  ^abcPabc  — .?  S  Pfjk  Pfan 

ljk  las 

ibcDai^bj^ck  +  (sl)abc^aiDbj  ^ck  +  (tx)abc^ai^bjDck] 


[(rx)abcDal^bai^cn  +  (Sx)abc^alDbm^cn  +  (ti)abc^al^bBDcn]  + 


(5.2.7) 


[(ry)abcDai^bj^ck  +  (sy)abc^ail^bj^ck  +  (ty)abc^ai^bjDck] 
[(ry)abcDai^bm<5cn  +  (Sy)abc^all^bB^cn  +  (ty)abc^al^>mDcn]  + 
[(rz)abcl-^ai^bj^ck  +  (Sz)abc^aiDbj^ck  +  (tz)abc^ai^bjl^ck] 
{(r|)abcDal^>ai^cn  +  (sljabc^alf^ba^cn  +  (tf)abc^al^bmD  “]]• 
Introducing  these  expressions  into  (5.2.2)  gives 
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J[P]' «  ZcM')Ms)M<-)  I  <iei[J«]  |  ,bc[-i  Sk  I/Sil  P!.„it 

f[(rx)abcDai^bi^ck  +  (Sxjabc^ai^bj  ^ck  +  (tx)abc^ai^bjDck] 
[(rx)abc^al^ba^cn  4  (sl)abc^alDbB^cn  4  (tx)abc^al^b«Dcn]  4 
{(ry)abcDai^bj5ck  4  (Sy)abc^ai^bj^ck  4  (ty)abc^ai^bjDck] 
[(ry)abcDal^ba^cn  4  (Sy)abc^al^bB^cn  4  (ty)abc^al<5baDcn]  4 
[(rl)abcDai^bj^ck  4-  (sl)abc^ai^bj^ck  4  (tt)abc^ai^bjDck] 
[(rz)abcDal^bm^cn  4  (s|)abc^a)Db«£cn  4  (tf)abc^al^b»Dcn]j 
4  V|bc-  ,?k^[(rx)abcDai^bj^ck  +  (Sx)abc^aiDbjfck  4  (tx)abcfci<VjDck]6x  4 
[(ry)abc^ai^bj^ck  4  (  S  y )  abc  <^a  b j  <5ck  4  (ty)abc^ai^bjDck]Cy  + 
[(rl)abcDai^bj^ck  4  (sljabc^ai^bj^ck  4  (tf)abc^ai<5bjDck]6zj  P  jjkj  drdsdt 

-  £  Pe.  (5.2.8) 

p  =  i 

Evaluating  the  integrals  by  introducing  wa=  f  "Va(r)dr  according  to  (1.5.16) 
gives 


J[P]e«-J  E  £  £  P?jk  PfBnAtwaWbWc|det[Je]|abc 
abc  ij  k  Inn 

|[(rx)abcDai^bj^ck  4  (Sx)abc^aiDbj ^ck  4  (tx)abc^ai^bjDck] 

[(rx)abcDal^bn^cn  4  (sl)abc^alDbB^cn  4  (tx)abc^aI^bniDcn]  + 

[(ry)abcDai^bj^ck  4  (Sy)abc^aiDbj^ck  4  (ty)abc^ai^bj^ck] 

[(ry)abcDal^bia^cn  4  (Sy)abc^alDbm^cn  4  (ty)abc^al^bB^cn]  4 

[(rz)abcDai<5bj  ^ck  4  (Sz)abc^aiDbj  $ck  4  (tz)abc^ai^bjDck] 

[(rz)abc^al^bB^cn  4  (Sz)abc^alDba^cn  4  (t|)abc^al^baDcn]j 

4  £  V|bc-  £  PijkWaWoWc|  det[Je]  |  abc 
abc  ljk 

^[(ri)abcDai^bj^ck  4  (Sx)abc^aiDbj  ^ck  4  (tl)abc^ai^bjDck]Cx  4 

[(fy)abcDai^bj^ck  4  (Sy)abc^aiDbj  ^ck  4  (ty)abc^ai^bj^ck]Cy  4 

[(rl)abcDai^bj^ck  4  (sl)abc^aiDbj ^ck  4  ( 1 1) abc ^a i ^bj D ck] j 

6 

-  £PC. 

p*i 


(5.2.9) 
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Excluding  the  surface  integral,  this  represents  the  governing  functional  in  terms 
of  the  local  co-ordinates. 

5.3  Surface  Integral 
The  surface  integral, 

E  P*  =  E  f 1  f 1  Pe(n-Ve)dAP,  (5.3.1) 

contains  contributions  from  each  of  the  six  sides  of  each  element.  We  must 
express  it  in  terms  of  the  unknown  expansion  coefficients  in  the  local  (r,s,t)  co¬ 
ordinate  system.  Rearranging  (5.3.1)  by  combining  the  surface  unit  normal 
vector,  n,  with  the  differential  area,  dA,  gives 

E  fiflpeye.dAP.  (5.3.2) 

p  —  1  p=K-K-I 

Introducing  these  expressions  for  each  of  the  six  faces  on  the  local  element  gives 

E  Pe  =  (5-3.3) 

p=i 

“/  'f 1  pe^e'(x3e*  +  yaey  +  zje^drds 
PeVe*(x2ex  +  y2 ey  +  Z2ez)drdt 

+J  jj*  |  P*Ve*(x3ex  +  y3ey  +  z^drds 

~f\f\  ^)e^re’(X2ex  +  z2ezMr^t 

-j  | J  j  peye-  (xtex  4-  yiey  +  z^dsdt 
+ /j/j  PeVe>(x1ex  +  yiey  +  z^dsdt, 
where  we  expressed  the  area  vectors  shown  in  (5.3.2)  in  terms  of  their 
component  parts  with  the  assistance  of  §A.7  in  appendix  A.  Introducing  the 
expansion  polynomials  into  these  surface  integrals  gives 

E  Pe  s  (5.3.4) 

p=i 

-  f 1  f !  E  PIbi^Ibi-(x3ex  +  y3ey  +  z3ez)abi^a(r)^(s)drds 

+  f 1  f 1  E  P|2cVI2c-(x2ex  +  y2ey  +  z2ea)a2c^a(r)^c(t)drdt 
J -\J  -  lac 
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+  f  lf 1  £  Pab2VIb2*(x3Cx  +  y&y  +  Z3ef)ab2lk(r)V>b(s)dids 
J -w  -i  ab 

~  J*1  f  1  £  PaicV|ic*(x2ex  +  y2®y  +  Z2CB)alc$a(r)V,c(t)dxdt 
-W  -1  ac 

-f  1f1  £  P1bc^!bc*(xiex  +  ytey  +  zie*)ibc^b(s)^c(t)dsdt 
»  - 1»  - 1  be 

+/1/ )  £  PSbc^1bc-(xiex  +  ytey  +  ziet)2bc^(s)^c(t)dsdt. 

»  -1*'  -1  be 

Evaluating  the  integrals  using  wa=/  ‘|^(r)dr  gives 

E  Pe  a  (5.3.5) 

p*i 

—  £  PlbiVfbr  (x3ex  +  ysCy  +  Z3e*)abiwawb 

ab 

+  £  PI2cVI2c-(xjex  +  y^y  +  z2ez)a2cwawc 

ac 

+  £  PIb2V|b2-(x3ex  +  y3ey  +  z3eJ.)ab2WaWb 

ab 

—  £  PlicVaic,(x2ex.+  y2ey  +  z2ez)aicwawc 
ac 

—  £  P!bcV!bc-(xiex  -f  y2ey  +  ziez)ibcwbwc 

be 

+  £  P5bcVSbc*(xlex  +  yiey  +  z1ez)2bcWbWc. 

DC 

w  w  v  v 

Finally,  performing  the  inner  product  using  Vlbc^UabcCx+VabcCy+WIbcez  gives 

EPe«  (5.3.6) 

p=i 

~  £  PIbi(Ulblx3  +  VIbiy3  +  WIbiz3)abiwawb 

ab 

+  £  Pa2c(U|2cX2  +  Va2cy2  4-  Wa2cz2)a2cWaWc 

ac 

+  £  Pab2<UIb2X3  +  V|b2ys  +  WIb2Z3)ab2wawb 
ab 

—  £  Paic(UaicX2  +  V|icy2  +  Waicz2)alcwawc 

ac 

—  £  P!bc(U?bcXl  +  Vlbcyi  +  W1bcZl)lbcWbWc 
be 

+  £  P5bc(USbcXi  +  VSbcyt  +  W5bcZi)2bcWbWc. 
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5.4  Variation  of  the  Functional 

All  of  the  pieces  corresponding  to  the  terms  in  the  governing  functional 
exist  in  terms  of  the  local  co-ordinates  and  expansion  functions  and  coefficients. 
Substituting  (5.3.6)  into  (5.2.9)  and  performing  the  remaining  inner  product 
gives 

J[P]e  8-i  S  E  £  P?jk  PS.nAtWaWbWc !  det[ Je]  |  abc 
abc  ij  k  Ian 

^[(rx)abc®ai^bj  ^ck  +  (si)abc^aiDbj  ^ek  +  (t|)abc^ai^bjDck] 
[(r|)abcDal^b*^cn  +  (Sx)abc^alDba^cn  +  (tx)abc^al^bmDcn]  + 
[(ry)abcDai^bj^ck  +  (sy)abc^aiDbj  ^ck  +  (ty)abc^ai^bjDck] 
[(ry)aboDal^ba^cn  +  (Sy)abc^all^ba^cn  +  (ty)abc^al^>»Dcu]  + 
[(rz)abcDai^bj^cV  +  (sf)abcfciDbj^ck  +  (tf)abc^ai^bjDck] 

[(rz)abcD  al  &bn  +  (Sz)abc^albbB^cn  +  (t|)abc^al^bnDcn]j 

+  JB.  P  ljkWaWbWc  |  det[Je]  |  abc  (5.4.1) 

|^[(rl)abcDai^bj  ^ck  +  (Sx)abc^ai^bj  ^ck  +  (t|)abc^ai^bjDck]U|bc+ 
[(ry)abcDai^bj  ft:k  +  (sy  )abc^aiDbj  ^ck  +  (ty)abc^ai^bjl^ck]V|bc+ 
[(rl)abcDai^bj  ^ck  +  (sl)abc^aiDbj^ck  +  (t|)abc^ai^bjDck]Wfbcj 

“  £  PIbl(U!biX3  +  V|biy3  +  W|biZ3)ablWaWb 

ab 

+  £  Pa2c(Ua2cX2  +  Va2cy2  +  Wf2cZ2)a2cva^c 

ac 

+  £  Pab2(Utb2*3  +  V|b2y3  +  W|b2Z3)ab2WaWb 
ab 

-  S  PIlc(UIicX2  +  Vllcy2  +  WI|Cz2)aicwawc 

ac 

-  £  P!bc(U!bcX|  +  V^bcyi  +  WtbcZi)tbcWbWc 

OC 

+  £  PSbc(USbcXi  +  VSbcyi  +  WSbCzi)2bcWbwc. 

DC 

The  extremum  of  J[P],  or  the  variation  5J[P]/5P,  equals 
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d3[?)*/d?  a  -  £  £  P?jkAtwaWbWc|det[Je]|abc 
abc  ij* 

^[(ri)abcDai^bj^ck  +  (Sx)abc^aiDbj^ck  +  (tl)abc^ai^bj'DckJ 

[(rI)abcCal^b«^cn  4"  (Sx)abc^al®b*^cii  4"  (tx)abc^al^)*l^cn]  4* 

[(ry)abcDai^bj  ^ck  4-  (Sy)abc^aiDbj  ^ck  4"  (ty)abc^ai^bj^ck] 

[(ry)abc^al^bB^cn  4"  (Sy)abc^al^bB^cn  4"  (ty)abcfcl^b«l^cn]  4* 

[(rl)abcDai^bj^ck  +  (sl)abc^aiDbj^ck  +  (tx)abc^ai^bjDck] 

[(rI)abcDal^b«<5cn  4"  (s|)abc^al^b*^cn  4"  (ti)abc^al^b«^cnl| 

+  £  waWbWc|det[Je]|abc  (5.4.2) 

abc 

|^[(ri)abcDai^bj  ^ck  4"  (Sx)abc^ai^bj^ck  4"  (tx)abc^i^bj^ck]Ufbc  4" 
[(ry)abcDai^bj  ^ck  4"  (Sy)abc^i^bj^ck  4"  (ty)abc^ai^bj'Dck]  Vabc  4" 
[(ri)abcDai^bj  ^ck  4"  (s|)abc^aiDbj^ck  4"  (t|)abc^ai^bj®ck]Wabcj 

-  £  (UlbiXa  4-  V|biy3  4-  W|biZ3)ablWaWb 
ab 

+  £  (Ua2cx2  +  Va2cy2  4"  W|2cz2)a2cwawc 

ac 

4-  £  (U|b2x3  4-  V|b2y3  4-  W|b2Z3)ab2WaWb 
ab 

-  £  (Uaicx2  4-  V|lcy2  4-  W|icZ2)aicWawc 

ac 

-  £  (UlbcXi  4-  Vtbcyi  4-  WfbczOibcWbWc 
be 

4-  £  (Ufbcxl  4"  V?bcyi  4"  W5bcZl)2bc^bwc- 
be 

The  middle  term  in  (5.4.2)  may  appear  more  conveniently  as 

£  wawbwc|det[je]|abc  (5.4.3) 

abc 

£[(*x)abcl^ai^bj^ck  4*  (sl)abcfcjl^bj^ck  4"  (t|)abc^ai^bjl^ck]Uabc  4" 
[(ly)abcDai^bj^ck  4-  (Sy)abc^aiDbj^ck  4-  (ty)abc^ai^bjDck]V|bc  4- 
((rl)abc-Dai^bj^ck  4-  ( S |) abc ^aiD bj ^ck  4-  (tf)abc^ai^bjDck]Wfbcj 
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=  S  V|bcijk-V|bcWaWbWc|det[Je]Ubc 

a  oc 

where  the  divergence  operator^  equals 

^abcijk  - 

[(rx)abc^ai^bj^ck  +  (sx)abc^aiDbj^ck  +  (tx)abc^ai^bjDek]ex  + 

[(ry)abcDai^bj^ck  +  (Sy)abc^ai^bj^ck  4*  (ty)abc^ai^bjHck]Cy  + 
[(rI)abcDai^bj  ^ck  +  (si)abcfci^bj^ck  +  (t|)abc^ai^>jDck]Cz- 
Introducing  the  quantity  nf.nabc>  defined  by 

Ilfanabc  =  — |  det[je]  ]  abcwawbWc^lafcb^nc> 

into  (5.4.3)  gives 


(5.4.4) 


(5.4.5) 

(5.4.6) 


^^Vtbcijk’ VabcWaWfoWc) det[Je]  |  abc - ^S^Vfunijk*  H5mnabcV|bo 

for  the  middle  term  in  (5.4.2).  The  free  set  of  indices  in  this  expression  equal 
ijk,  while  those  in  the  first  term  in  (5.4.2)  equal  lmn.  Rearranging  the  indices 
gives  an  equivalent  form  for  the  right-hand  side  of  (5.4.6),  written  as 


,?k  a|jc  ^Tjklan '  HijkabcVfbc- 


(5.4.7) 


The  extremum  of  J[P]  occurs  when  the  expression  in  (5.4.2)  equals  zero. 
Substituting  (5.4.7)  for  the  middle  term  and,  setting  <5J[P]/<5P=0  gives 


AfanijkPijkM  —  — Vijklan  *  H' ljkabcVfbc  +  Se, 


ijk  ’  •  ij 

where  the  surface  integral  appears  as 


Se  = 


(5.4.8) 

(5.4.9) 


W  v  w 

~  E  (UIbix3  +  V|biy3  +  W|biz3)abiwaWb 

ab 

+  E  (Uf2cx2  +  V|2cy2  +  WI2cz2)a2cwawc 

ac 


JThe  term  Vabcijk*  Vabc,  obtained  from  the  expression  V-V  does  not  equal  the 
divergence  of  velocity.  We  write  the  divergence  of  velocity  V-V  as  VabcijkVijk 
Hence,  the  order  of  indices  is  important. 
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+  £  (U|b2X3  +  VIb2y3  +  WIb2Z3)ab2WaWb 

ab 

-  £  (U£icX2  4-  V|icy2  4-  WIicz2)aicwawc 
ac 

-  £  (U£bcxi  4-  Vfhcyi  +  ^tbcZi)ibcWbWc 
be 

4-  £  (USbcXi  +  V^hcyi  +  W1bczi)2bcwbwc> 
be 

and  the  quantity  Atmnijk  as 

A!«nijk  =  Jc  wawbwc|  det[Je]  I  abc[  (5.4.10) 

[(rx)abcDai^bj^ck  4"  (Sx)abc^ai^bj^ck  4"  (tx)abc^ai^bjBck] 
[(rf)abcDal£b»^cn  4-  (si)abc^alDb»^cn  4*  (t|)abc^al^b»Dcn]  4* 
[(ry)abcDai^bj^ck  4-  (Sy)abc^ai^bj^ck  4"  (ty)abc^ai^bj^ck] 
[(ry)abcDal^bm^cn  4"  (Sy)abc^alBbn^cn  4"  (ty)abc^al^bB^cn]  4* 
[(ri)abc^ai^bi^ck  4"  (si)abc^ai^bj  ^ck  4”  (tl)abc^ai^bj^ck] 
[(rI)abcDaJ^ba^cn  4"  (sf )abc^alDba Sea  4-  (tf)abc^al^bsiDcn]j  • 


5.5  Solution  for  Pressure 

Equation  (5.4.8)  applies  to  a  single  element.  Adding  the  contributions 

from  all  elements  gives  (5.5.1) 

E  E  t  E 

£  £  Afnnijk-P  ijkAt  =  £  £  £  V  jjklmn  *  HijkabcVabc  4-  £  Se. 

e  =  l  ijk  e*l  ijk  abc  e=l 

In  this  expression,  the  pressure  coefficients,  Pfjk,  remain  unknown;  the  velocity 


coefficients,  V|bc,  known  from  the  non-linear  step;  and  the  surface  integral 

unknown,  since  it  contains  unknown  velocity  coefficients  V|bc-  By  summing 
over  all  elements,  and  taking  advantage  of  the  continuity  of  the  velocity  on  the 
faces  shared  between  elements,  we  greatly  simplify  the  expression  for  the  surface 
integral.  Since  the  surface  integrals  contain  velocities  and  not  derivatives  of 
velocity,  and,  since  the  velocity  remains  continuous  throughout  the  domain,  the 
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terms  cancel  on  faces  shared  between  two  elements.  Thus  the  summation  over  all 
elements  yields  a  contribution  from  only  those  elements  on  the  boundary  of  the 
domain.  Unfortunately,  the  surface  integral  contains  unknowns  which  represent 
the  velocity  after  the  pressure  step.  Hence,  (5.5.1)  contains  two  unknowns: 

pressure  P!jk  and  velocity  V?jk-  Along  surfaces  where  we  specify  essential 
boundary  conditions  on  velocity,  the  surface  integral  depends  on  known 
quantities;  along  surfaces  where  we  specify  natural  boundary  conditions,  the 

velocity,  V?jk,  remains  unknown.  We  may  circumvent  this  problem  by 

approximating  the  surface  integral  of  n*V  along  natural  boundaries  by  n*V, 

since  we  know  V  from  the  non-linear  step.  However,  since  V  does  not  satisfy 

boundary  conditions,  a  significant  error  may  accrue  from  these  integrals; 

furthermore,  these  integrals  over  the  boundaries  wher^  both  natural  and  essential 

conditions  occur  may  require  extensive  computational  time.  Therefore,  we  use  an 

entirely  different  technique  for  the  surface  integral  expression. 

E  w 

Since  E  Se  represents  the  net  flux  of  velocity  leaving  the  domain, 
e  =  l 

/Wtt'V)dA,  we  choose  to  compute  the  integral  of  n*V  over  the  inlet  and 
exit  boundaries  and  apply  the  difference  of  these  terms  to  the  exit  velocity  in  the 
form  of  a  correction.  Hence,  no  longer  do  we  satisfy  (5.5.1)  exactly  at  each  of 
the  nodal  points.  Instead,  (5.5.1)  holds  only  in  a ’'global"  sense.  We 

approximate  the  surface  integral  by 

E  -  _ 

S  S*  «  Vout'nAoat 

<5-5-2) 

where  V0ut  equals  the  average  outflow  velocity  correction  and  Aout  equals  the 
outflow  area.  The  velocity  correction  applies  to  the  velocity  after  the  non-linear 
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step  according  to 


(Vijk)corrected  —  Vtjk  —  Vout>  (5-5.3) 

where  node  ijk  and  element  e  correspond  to  those  nodes  looted  along  the  exit 

portion  of  the  domain,  Vtjk  represents  the  known  velocity  after  the  non-linear 
step,  and  Vout  the  average  velocity  from  (5.5.2).  We  use  this  corrected 
velocity  in  (5.5.1)  in  place  of  the  original  uncorrected  velocity,  and  we  drop  the 

surface  integral.  The  final  form  equals  (5  5  4) 

e  e  - 

5j  S  AfanijkPijk^t  =  E  S  S  Vtjkl*n*ntjkabc(Vabc)corrected- 
e  =  l  ijk  J  e=l  ijk  abc 

Thus,  we  achieve  the  standard,  linear  form,  Ax=B,  for  the  unknown  pressure 
coefficients.  After  computing  the  pressure,  the  velocity  after  the  pressure  step 
results  from 

V?jk  =  %-njkPSjkAt,  (5-5.5) 

equivalent  to  the  discritized  form  of  (3.2.2)  applied  at  node  ijk  of  element  e. 

Incidentally,  a  proviso  must  accompany  quantities  obtained  from 
expressions  with  derivatives,  such  as  (5.5.5).  Since  only  the  variables  appearing 
in  the  essential  boundary  conditions  remain  continuous  across  elements,  and  not 
their  derivatives,  the  gradient  of  pressure  contains  different  values  at  the 
intersection  between  two  elements-one  value  from  each  of  the  two  elements. 
Since  only  one  value  may  exist  at  any  point  in  the  global  numbering  scheme, 
before  equations  such  as  (5.5.5)  get  solved,  we  must  combine  the  two  values  for 
VP  at  node  ijk  by  taking  the  average  value  from  the  two  elements.  This  seems 

innocuous  enough;  unfortunately,  we  break  the  zero  divergence  rule  on  Vtjk  in 
the  process.  This  may  induce  important  errors  in  certain  problems. 


CHAPTER  VI 
VISCOUS  STEP 


6.1  Variational  Approach 

The  governing  differential  relation  for  the  viscous  step  equals  Helmholtz’s 
equation,  given  by 

inO.  (6.1.1) 

We  employ  techniques  from  calculus  of  variations  to  generate  a  matrix  expression 
for  the  unknown  expansion  coefficients  $yk.  The  goal  remains  finding  a 
functional  yielding  Helmholtz’s  equation  after  application  of  the  Euler-Lagrange 
formula.  Assume  for  the  moment,  that  the  functional  whose  Euler-Lagrange 
equation  corresponds  to  (6.1.1)  with  homogeneous  boundary  conditions  given  by 

o4  +  0(n-V)$  =  O,  on  dft,  (6.1.3) 

equals 

J»]  =  Jfl  [-igrad*:giad*-  JAJ+-t-+.F]dV.  (6.1.2) 

This  functional  applies  to  domains  with  rigid  boundaries.  In  fluid  dynamics, 
both  homogeneous  and  non-homogeneous  boundary  conditions  arise.  In  general, 
the  two  possible  types  of  non-homogeneous  boundary  conditions  equal  Dirichlet 
(essential)  when  /3=0,  given  by 

*=a,  (6-1.4) 

and  Neumann  (natural)  when  a=0,  given  by 

(n.V)<fr=6-  (6.1.5) 

We  may  incorporate  these  non-homogeneous  boundary  conditions  in  the  standard 
functional  by  adding  a  boundary  integral. 

Since  $  must  remain  fixed  along  surfaces  containing  non-homogeneous 
Dirichlet  boundary  conditions  according  to  (6.1.4),  the  variation  in  written 
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as  $*=4+e^,  must  equal  4-  other  words,  we  impose 

tfdfl)  =  0  (6.1.6) 

on  surfaces  with  Dirichlet  boundary  conditions.  We  impose  no  restrictions  on  j>, 
along  boundaries  with  natural  conditions.  Instead,  the  natural  or  free  boundary, 
where  4  may  vary,  requires  the  standard  constraint,  n-dF/d(grad$)=0,  where 
F  equals  the  integrand  of  the  functional.  The  result,  after  applying  this 
condition  to  the  functional  given  in  (61.2),  equals  n«grad$=0,  which  does  not 
satisfy  the  non-homogeneous  natural  boundary  conditions.  We  overcome  this 
problem  by  adding  a  surface  integral  to  the  functional  giving 

=  Jq  [-  *grad$:grad$  -  -  $-F]dV  + 

JdCln  (6.1.7) 

where  dQn  represents  the  surface  containing  natural  boundary  conditions.  The 

variation  of  the  surface  integral  contributes  a  term,  which  when  combined  with 
n-3F/5(grad$),  yields  the  proper  non-homogeneous  natural  boundary  conditions. 
Until  now,  we  accepted  the  claim  that  the  functional  shown  in  (6.1.7)  yields 
Helmholtz’s  equation  with  appropriate  boundary  conditions  without  proof.  We 
must  show  that  the  extremum  of  this  functional  does  indeed  yield  Helmholtz’s 
equation  and  boundary  conditions  for  the  viscous  step. 

We  begin  by  considering  a  functional  represented  as 

■>[♦]  =  Jn  n*.  ♦.  grai«dV.  (6.1.8) 

The  variation  in  this  functional  equals  the  difference  between  the  functional 
evaluated  at  and  $.  Thus,  &J=J[<fr*]-J[$]  may  appear  as 

=  Jq  [F(x,  <J>  +  grad$  +  cgrad^) 

-  F(x,  $,  grad$]dV.  (6.1.9) 

Expanding  the  integrand  in  a  Taylor  series  gives 
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63  =  Jfl  +  'flg?ad$c:grad^dV’  (6.1.10) 

which,  after  application  of  the  chain  rule,  yields  > 


Applying  this  expression  to  the  functional  in  (6.1.7)  results  in 

a  =  efQ  [V  -grad*  -  A2*  -  F]  -  *1V  +  efQ  V  •  [-  grad$.  tfdV 

+  cJ^flng*lW(T.  (6.1.12) 

The  divergence  theorem  allows  us  to  express  the  second  integral  on  the  right  as 
cJfi  =  ~  efd n  n*(grad^*i&)da.  (6.1.13) 

The  integrand  on  the  right-hand  side  equals  zero  along  the  boundary  where 
essential  conditions  apply,  since  ^=0  in  those  regions;  the  only  contribution 
comes  from  boundaries  containing  natural  boundary  conditions.  Using  the 
natural  boundary  conditions  given  in  (6.1.5)  gives 

f  Jp  V-I-  grad$*$dV  =  -  c g-i^a  (6.1.14) 

which  exactly  cancels  the  last  term  in  (6.1.12).  Thus,  the  variation  in  J  equals 

fiJ  =  €jfj[V-grad^-A^-F]-^dV.  (6.1.15) 

At  an  extremum,  we  require  &J=0  for  all  admissible  rfr,  in  particular,  we 
require  &J=0  for  all  admissible  ^  which  vanish  on  the  surface  containing 
essential  boundary  conditions.  Because  of  the  arbitrariness  of  $  inside  fl, 

&J=0  implies 

V2<^-A^-F  =  0,  in  Q.  (6.1.16) 

This  equals  Helmholtz’s  equation,  confirming  our  supposition. 

As  indicated  previously,  the  variational  form  incorporates  the  natural 
boundary  conditions  in  the  functional,  while  the  essential  boundary  conditions  do 
not  appear.  The  functional  expression  makes  no  reference  to  a  particular  domain; 
hence,  it  applies  to  both  a  single  element  or  the  entire  domain.  In  the  following 


sections  we  express  this  functional  in  terms  of  the  local  co-ordinates  on  a  single 
element.  Each  element  contributes  both  a  volume  integral  and  a  corresponding 
surface  integral.  A  functional  applicable  to  the  entire  domain  results  from 
summing  over  all  elements. 

6.2  Representation  in  Local  Co-Ordinates 

The  governing  functional  written  in  global,  or  physical  space,  equals 

J[«  =  X  [-isra#gra<4-iAiH-fF]dV 

+  Jan”  (6-21) 

We  must  transform  the  functional  expression  from  the  global  (x,y,z)  co-ordinate 
system  to  a  local  (r,s,t)  system,  where  elements  appear  as  cubes  whose  local  co¬ 
ordinates  range  between  ±1.  The  functional  representation  for  a  single  element  in 
local  co-ordinates  equals 

J[<j»]e  =  / 1/  j/ 1(-  igrad$:grad$  -  $ A2^-^  -  $-F]e  J  det[Je]  |  drdsdt 

(Consult  §A.3  for  details  of  the  differential  volume  conversion  between  global 
and  local  co-ordinates.)  The  surface  integral  expresses  a  summation  over  the  six 
element  faces.  Remember,  this  term  applies  only  along  surfaces  characterized  by 
natural  boundary  conditions.  Notice  that  the  local  co-ordinate  system  occupies  a 
position  within  an  element  such  that  only  two  co-ordinates  vary  along  a  given 
face,  while  the  third  co-ordinate  remains  fixed  at  ±1.  The  gradient  operator  in 
the  global  system  equals 

grad$e  =  Hex  +  |yey  +  !!  ^  (6.2.3) 

In  the  local  co-ordinate  system,  [r(x,y,z),  s(x,y,z),  t(x,y,z)],  this  operator  equals 

grad*'  =  [ffr!  +  f'sS  +  tS]e*  (6.2.4) 
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+  lg  rf  +  |f»?  +  |f  t?)**  +  I#'*?  +  #'*?  +  |f«3*> 

where  the  partial  derivatives  of  the  local  co-ordinates  with  respect  to  the  global 
ones  appear  according  to  (A.  2.7)  as 

_  /YsZt"  JtZs\e.  re  _  /*tzs“  xsZne.  Te  _  f*s Yt~  *tJse. 

~  ^  clit[JeH  ’  y  "  (  ditpeH  ’  2  ~  (^etpfer  : 


r| 


*? = ( Wf )e;  s' - W1  ■*  -  ( WP)e; 

*?  =  P&cjtfr-  *?  = 


(6.2.5) 


(6.2.6) 


(6.2.7) 


The  Jacobian  of  the  transformation  equals 


Je  = 


f  Xr  Xs 


Xt 


yr  ys  yt 

zr  Zs  ZtJ 


(6.2.8) 


according  to  §A.2.  Subscripts  indicate  differentiation  with  respect  to  local  co¬ 
ordinates;  superscripts  define  the  particular  dement  under  study. 

6.3  Approximating  Functions  in  Spectral  Space 

The  discrete  expansion  of  a  function  u(r,s,t)  in  terms  of  a  finite  sequence 
of  orthogonal  functions  of  n-th  degree  approximately  equals 


u(r,s,t)e  a  £  £  £  u!jkV'i(r)V>j(s)V’k(t),  (6.3.1) 

irQ  J  5 0  1*0 

where  the  scalars  ufjk  equal  the  discrete  expansion  coefficients  of  u(r,s,t) 
rdative  to  the  basis  Notice  that  this  expression  is  not  of  n-th  degree 

unless  l=m=n.  As  indicated  in  §1.5,  differentiation  may  take  place  in  either 
spectral  space  or  physical  space.  Normally,  we  use  differentiation  in  physical 
space.  The  derivative  with  respect  to  the  local  co-ordinate  r  in  the  local  space 
equals 


Ur{r,S,t)e 


t?o  j?0  kf0^jk|^i(r)^j(s)^k(t). 


(6.3.2) 


52 


We  do  not  compute  derivatives  in  this  manner  since  this  would  entail  computing 
analytical  expressions  for  the  derivatives  of  the  basis  functions;  instead,  we 
compute  derivatives  only  at  the  nodal  points  (ra,Sb,tc).  This  gives 

Ur(ra}Sb,tc)e  *< 


E  E  E  u?jkMra)^j(»b)^(tc).  (6.3.3) 

1  =  0  j=0  k  =  0 

Since  the  basis  functions  satisfy  both  V'i(ra)=  (5ai  and  Dai=|&(ra),  this 
expression  reduces  to 


ur(ra,Sb,tc)e  »  ,SQ  ,SQ  ^ufjkDai^j^k,  (6.3.4) 

which  may  be  written  more  compactly  as 

(Ur)abc  *  ,?k  fiijk  Dai^bj^ck*  (6.3.5) 

Similarly,  derivatives  with  respect  to  the  s  and  t  directions  equal 

(us) abc  *  Sk  Uijk  <5aiDbj^ck  (6.3.6) 

and 

(ut) abc  ~  ,|k  fifjk  ^ai^bjDck,  (6.3.7) 

respectively.  The  function  u(r,s,t)  approximated  at  the  node,  (ra,Sb,tc)  equals 


u(ra,Sb,tc)e  « 

_E  E  E  Uijk^,i(ra)^'j(sb)^,k(fc)  (6.3.8) 

1  =  0  j  =o  k  =  0 

or,  in  compact  form, 

Uabc  ~  .?  Uijk^ai^bj^ck-  (6.3.9) 

ljk 

With  these  expressions  for  functions  and  their  derivatives,  all  the  ingredients 
necessary  to  represent  the  functional  in  variational  form  exist.  Incidentally,  the 
first-order  variational  form,  as  opposed  to  the  second-order  differential  form,  does 
not  require  a  second  order  derivative,  clearly  an  advantage  given  the  complexity 
of  the  first  order  derivatives. 

6.4  Approximation  of  Products  of  Functions 
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Since  the  variational  form  of  our  functional  contains  products  of  functions, 
some  of  which  are  known  (e.g.,  the  Jacobian  or  the  body  force),  while  others, 
such  as  velocity,  remain  unknown,  we  need  a  system  to  express  products  of 
functions.  The  two  most  common  methods  of  expressing  the  product  of  functions 
u  and  v  equal  the  product  of  the  series, 

u-v  a  (Inu)(I„v),  (6.4.1) 

or  the  series  of  the  product, 

u-v  »  In(u«v).  (6.4.2) 

The  first  form  yields  a  polynomial  with  nJ  terms  since  it  involves  term-by-term 
multiplication  of  two  series;  the  second  method  involves  a  polynomial  with  n 
terms  since  the  product  u  and  v  occurs  before  the  expansion  in  terms  of  the 
nodal  coefficients.  We  use  the  second  method.  Using  the  compact  notation,  the 
product  of  u  and  v  equals 

ue-ve  »  Ek  u1jk-v?jk^i(r)V»j(s)t^(t).  (6.4.3) 

A  natural  extension  of  this  expression  allows  a  representation  of  products  of  three 
or  more  terms. 


6.5  Surface  Integrals 

The  term  representing  the  surface  integral  in  (6.2.2)  equals 
Ie  =  EIJ*=  E  pp$*-g*dAP. 

p  =  i  p  =  i-'  - 1-/  - 1 

Using  the  inhomogeneous  boundary  condition,  g=(n-V)$,  gives 


which  after  rearranging  gives 


Ie=  lJ  j  <|>e(nedAp-V)$e. 


(6.5.1) 


(6.5.2) 


(6.5.3) 


The  projection  of  the  differential  surface  area  in  the  direction  of  the  outward  unit 
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normal  on  each  of  the  six  faces  equals 

C*  Cy  ©2  e 

nedA*  =  —  xryrZrdnis,  (6.5.4) 

.  xs  ys  z8. 

at  t=— 1; 

Cx  Gy  ez  c 

nedA2  =  xr  yr  zr  drdt,  (6.5.5) 

.  yt  zj 

at  s=-l; 

Ci  6y  Cz  e 

nedA3  =  x  yr  zr  drds,  (6.5.6) 

.  xs  ys  zs. 

at  t— 1; 

Cx  Cy  Cz  e 

nedA4  —  —  xr  yr  zr  drdt,  (6.5.7) 

.xtytzt. 

at  s=l; 

Cx  Cy  ^  e 

nedA5  =  -  xs  y$  zg  dsdt,  (6.5.8) 

.  xt  yt  zj 

at  r=-l;  and 

Cx  Cy  ez  e 

nedA6  =  xs  ys  zs  dsdt,  (6.5.9) 

.  xt  yt  zt. 

at  r=l.  (See  §A.7  for  details.)  Taking  the  inner  product  of  the  projection  of 
the  differential  area  in  the  direction  of  the  outward  normal  with  the  gradient 


operator  for  each  of  the  six  faces  gives 

n'dA'-V'  =  -  (»l|j  +  ‘yjy  +  t4)det[Je]drds,  (6.5.10) 

on  face  1; 

nedA2*Ve  =  —  (slgj  +  Sy^y  +  s||^)det[Je]drdt,  (6.5.11) 

on  face  2; 

nedA3*Ve  =  (tf^  +  ty^-  4-  t|^)det[Je]drds,  (6.5.12) 

on  face  3; 

nedA4*Ve  =  ^  +  sl^)det[Je]drdt,  (6.5.13) 
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on  face  4; 

n'dA >-r  =  -  (r||j  +  +  t$^)det[J']dsdt,  (6.5.14) 

on  face  5;  and 

n'dA‘.r=  (tl^  +  +  r^)det[J')dsdt,  (6.5.15) 

on  face  6,  where  the  derivation  of  such  quantities  as  Xx=dt/dx  appear  in  §A.2. 
Substituting  these  expressions  into  (6.5.3)  gives  for  each  of  the  six  faces: 

Ile  =  - ♦'(‘Sjf*  +  +  »8*V‘  [j'ldrds, 


on  face  1; 


on  face  2; 


on  face  3; 


on  face  4; 


on  face  5;  and 


on  face  6. 


iM = ♦'(«?§'  +  + ss§e)d<*mdrdt, 

i!e  -  + t5#')det[j']drds, 

I‘e=  +  sl|')det[Je]drdt, 

I*  =  ♦ew|'  +  rf§'  +  r|||e)det(J']<isdt, 

1 

I"  =  /.'/I  +'(r^‘  +  'yf'  +  r?§e)det(J')dsdt, 


(6.5.16) 


(6.5.17) 


(6.5.18) 


(6.5.19) 


(6.5.20) 


(6.5.21) 


Substituting  the  expression  for  the  derivatives  and  ||  from  (6.2.4) 

into  the  surface  integral  for  face  1,  gives 

1“  =  If'  +  stf  +  ‘I  £f>  +  ‘K'y  |f'  +  «• 


W  +  *» 


+  t|(rl  ff  +  s|  H'  +  t?  |f)fydet[J']drds. 


(6.5.22) 


We  must  introduce  the  expansion  polynomials  and  coefficients  into  these  integral 
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expressions.  The  result  for  face  1  equals 


(6.5.23) 


+  t|(r|^  +  fl '^)]abidet[Je]abI^ablV,a(r)V,b(s)<ir<is> 


assuming  the  quantity  within  square  brackets  refers  to  the  nodal  points  at 
ra,Sb,ti.  Similar  expressions  apply  for  faces  2  through  6.  Substituting 
expressions  for  the  partial  derivatives  of  $  from  (6.3.5),  (6.3.6),  and  (6.3.7) 

^ves  (6.5.24) 

Ile  ~  —f  1  f  1  E  .?  [t|(rxDai5bj^lk  +  Sx^ai^bj^lk  +  t|6ai^bjDlk)  + 

*  -K-l  ab  ijk 

ty(ryDai<$bj^lk  "b  Sy^aiDbj^lk  +  ty(5ai<5bjD  ik)  +  tf(r|Dai4.j  <5ik  + 

sf^aiDbj^ik  +  t|5ai^bjDik)]abi$iik$Ibidet[Je]ab#a(r)V,b(s)drds. 

The  only  variables  remaining  within  the  integrand  equal  V*(r).  i^s),  and  ip(t)\ 
therefore,  integration  reduces  to  /ii$i(r)dr=Wi.  Substituting  these  expressions 
into  (6.5.24)  gives  (6.5.25) 

IlB  “  "fb  ifJ^rSDaiMik  +  sf&iDbj*ik  +  ti<5ai^bjDlk)  + 
ty(ryDai^bj^ik  +  Sy^aiUbj ^lk  ■+■  ty^ai^bjD ik)  +  t|(r|Dai4)j  ^lk  + 

sf^aiDbj^ik  +  fz^ai^bjDik)]abl^fjk$abldet[Je]ablwawb- 
The  superscript  tilde  identifies  unknown  discrete  expansion  coefficients.  The 
other  terms  do  not  contain  the  tilde  (e.g.,  t|,  rl  and  det[Je])  since  they 
represent  known  continuous  functions  which  do  not  require  expression  in  terms  of 
approximate  expansion  coeffidents-these  coeffitients  depend  on  evaluation  of 
functions  at  the  nodes.  For  completeness,  the  surface  integrals  for  faces  2 
through  6  equal 

I2e  a 

—  E  E  [si(r|Dai62j£ck  +  s|5aiD2j^ck  +  txfci&jDck) 

ac  ijk 

sy(ry^*ai^2j^ck  +  Sy6aiD2j£ck  ty^ai^2jDck)  +  s|(r|Dai^2j^ck  + 


(6.5.26) 


sf£aiD2j<?ck  +  t|^ai^2jDck)]abl4ijk$i2cd6t[Je]a.2cWaWc; 

I3e« 

.?k[t!(r!Dai5bj«2k  +  Sx^aiDbj^2k  +  tI5ai^bjD2k)  + 
tyfry^ai^bj^k  +  Sy^aiDbj $2k  4*  ty£ai^bjD2k)  +  t|(r|Dai^bj52k  4" 
Sz$aiDbj^2k  4-  t  f  £a  i  5bj  D  2k)]  ab  1$  1  j  k$ab  2det  [  Je]  ab  2waw  b  > 

I4e  a 

E^  EJsKrSDai^ck  4*  sl^aiDjj^ck  +  tx^ai^ljDck)  4- 
sy(ry®ai^lj^ck  4-  Sy^aiDij^ck  +  ty£ai^ijDck)  4-  s|(r|Dai5ij  6ck  + 
sf^aiDij^ck  4"  t|5ai^ljDck)]abt$ijk^licdet[Je]aicwawci 

I5e« 

~  be  jj^tr*(r*^li^bi^ck  Sx^liDbj^ck  4-  tl^ii^bj^ck)  4" 
ry(ry^ii^bj^ck  4"  Sy^jjDbj^ck  4-  ty^H^bjDck)  4-  r|(r|Dii5bj5Ck  4- 
si^iiDbj^ck  4-  t|^ii5bjDck)]ibc^fjk^1bcdet[Je]ibcWbWc',  and 

I6e  ~ 

E_  .?k[ri(r*I>2i5bj ^ck  4-  S|^2iDbj^ck  4"  txfoi^bjDck)  4" 
ry(ryD2i^bj <5ck  4-  Sy^jDbj ^ck  4-  ty^i^bjDck)  4-  r£(r|D2i^bj5ck  4- 
s|^2il^bj^ck  4"  ti^2i^bjDck)]2bc$?jk$5bcd6t[Je]2bcwbWc- 
6.6  Application  to  Variational  Form 

The  governing  functional  approximated  in  discrete  form  equals 

■J  [4*3 6  tt  f  \f\J  ;ay~  ^abc^ibc:^abc^abc 

~  i^^Ibc^abc  ~  <^abc‘  Fabc]V,a(r)^fb(s)V,c(t) |  det[Je]  j  abcdrdsdt 

4-  E  Pe. 
p  =  i 

The  first  term  within  the  integral  on  the  right-hand  side  appears  as 

^abc^abc^abc^abc  = 

[(^)abe(rx)Ibc  4-  {^)abc(Sx)abc  4-  (j^)abc(tx)tbc]£x 


(6.5.27) 

(6.5.28) 

(6.5.29) 

(6.5.30) 


(6.6.1) 

(6.6.2) 
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^)abc(ry)fbc  +  (Jg)abc(Sy)abc  +  (^)abc(ty)abc]Cy 
^)abc(rz)abc  +  (^)abc(Sz)Ibc  +  (^)abc(tz)abc]Czj  * 
^)abc(rx)abc  +  (^)abc(Sx)abc  +  (^)abc(tx)fbc]Cx 
^)abc(ry)lbc  +  (^|)abc(Sy)abc  4-  (^)abc(ty)!bc]Cy 
ij)abc(rz)abc  +  (^)abc(&z)lbc  4-  (^)abc(tz)abc]Gzj  • 


Performing  the  inner  product  and  introducing  the  discrete  expansions  for  the 
partial  derivatives  gives 


Vfbc$ibc:^abc$abc  — 


(6.6.3) 


)abcDai^bj  ^ck  4*  (Sx)abc^aiDbj  ^ck  4*  (tx)Ibc$ai^bjDck] 
[(rx)IbcDal^bm<5cn  +  (Sx)abc^all^ba^cn  +  (tx)lbc^al^bmDcn]  + 
[(ry)abcDai^bj^ck  4*  (Sy)Ihc^ai^bj  ^ck  +  (ty)fbc^ai^bjDck] 
[(ry)abcDal^ba^cn  +  (Sy)Ibc^alDba^cn  4*  (ty)abc$al^b«Dcu]  + 
[(rz)abcDai^bj^ck  4-  (Sz)abc^aiDbj  ^ck  4"  (tz)abc^ai^bjDck] 
l(rz)lbc^al^ba^cn  4"  (Sz)abc^al^bn^cn  4-  (tz)Ibc^al^bmD  cn]]- 
Substituting  this  relation  into  (6.6.1)  gives 

w  *  ^woiw^oitepiikbcj-  *sk  ^k-*?.. 

|^[(rz)lbcDai^bj  ^ck  +  (Sz)abc^ai^bj^ck  +  (tx)abc^ai^bjDck] 
[(rx)abcDal^bis^cn  4*  (Sx)Ibc^alDb*^cn  4-  (tx)abc^al^b*f^cnj  4- 
[(ry)abcDai^bj^ck  4"  (Sy)lbc^ai^bj^ck  4*  (ty)abc^ai^bjDck] 
[(*y)abcDal^ba^cn  4"  (Sy)abc^alDbm^cn  4-  (ty)fbc^al^bi»Dcn]  4- 
[(rz)abcDai^bj^ck  4-  (Sz)lbc^aiDbj^ck  4-  (tz)Ibc^ai^bjDck] 
[(rz)abcDal^bn^cn  4-  (Sz)abc^alDbn^cn  4-  (tz)abc^al^bBiDcn]| 

—  -JA^Ibc'^abc  —  $Ibc' Ffbcl  drdsdt  4*  E  I^e.  (6.6.4) 

J  p  =  l 

Again,  the  only  variables  within  the  integrand  equal  $ r),  ip( s),  and  ip(i). 
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Using  (1.5.16)  for  /-i^i(x)dx  gives 

J[$]e  *  -  ,?k  -$!«n  I  det[Je]  |  abcWaWbwc 

^[(rx)IbcUai^bj  ^ck  +  (Sx)abc^aiUbj^ck  +  (tx)abc^ai^bj^ck] 
((rx)abcUal^b»^cn  +  (Sx)Ibc^alUbB^cn  +  (tx)abc^al^b«Ucn]  + 
Kry)IbcUai^bj^ck  +  (Sy)abc^aiDbj^ck  +  (ty)Ibc^ai^bjUck] 
[(ry)IbcUal^b«^cn  +  (Sy)abc^alUbB^cn  +  (ty)abc^al^b»Ucn]  + 
[(rz)abcDai^bj^ck  +  (Sz)abo^aiUbj^ck  +  (iz)abc^ai^bjUck] 
[(rz)abcUal^bm^cn  +  (Sz)abc^alUbo^cn  +  (tz)Ibc^al^b«iUcn]j 

-  $A2  £  ^.n^fBn  |  det[Jel  |  l*nWlWBWn 

inn 

-  s  ^fBn-FfMldetlJ^llBmWlW.Wn  +  £  Ipe. 

Imn  p=l 


(6.6.5) 


The  equation  which  minimizes  the  functional  J[$]  occurs  when  6J/fy\ma  =  0, 


and  equals 


l$lnm  ~  -  abc  i j  k  ^jk * det^  ® abcWaWbWc  (6-6-6) 

jji(rx)abcDai£bj5ck  +  (Sx)Ibc^aiUbj  ^ck  +  (tx)abc^ai^bjUck] 

[(rx)abcDal^b«^cn  +  (Sx)abc^alUbm^cn  +  (tx)ibc^al^bmUcn]  + 

[(ry)abcUai^bj  ^ck  +  (Sy)lbc^aiUbj^ck  +  (ty)abc^ai^bjUck] 

[(ry)abcUaI^bm^cn  +  (Sy)tbc^alDba^cn  +  (ty)abc^al^bnUcn]  + 

[(rz)lbcUai5bj  £ck  +  (Sz)lbc^aiUbj^ck  +  (tz)abc^ai^bjUck] 

[(tz)lbcUal^bm^cn  +  (Sz)abc^alUbm^cn  +  (fz)abc^al^bm^cn]j 

“  A2$?Bn  |  det^]  |  lnnWlWBWn  -  FfBn  |  detfJ6]  |  lBnWlWBWn  + 

SfdQ  n'^gfla$dcr  + f$(  p?i pe]  • 

The  last  two  terms  cancel.}:  Setting  fiJ/%Bn=0  gives  a  more  compact  form  of 


JThe  integrand  of  the  second-last  term  equals  n*5F/<5grad$  which  also  equals 
- n-grad$  for  the  functional  defined  in  (6.2.2).  The  integration  of  this  term 
takes  place  over  both  the  essential  and  natural  boundaries.  The  variation  of  the 
surface  integral  with  respect  to  $  equals  g  integrated  over  the  boundary  where 
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(6.6.6): 


.?^[A.fmnijk  —  ^2H?smijk]$ijk  —  .?^nf*njjkF^jk,  (6.6.7) 

where  II  and  A  appear  as 

nfnnijk  =  ~|  det[Je]  |  ijkWiWjWk^ij^j  5nk  (6.6.8) 

and  Af.nijk  s  JJ  det[Je]  |  abcwawbwcj  (6.6.9) 

[(rx)abc^ai^bj^ck  +  (Sx)&bc£aiDbj^ck  +  (tx)abc^ai^bj^ck] 

[(rx)ibcDal<5bB<5cn  +  (Sx)Ibc^alDbB^cn  +  (tx)abc^al^b*Dcn]  + 

[(ry)abcDai5bj  <5ck  +  (Sy)abc^aiDbj^ck  +  (ty)abc^ai^bjHck] 

[(ry)abcDal^bm^cn  +  (SyJIbc^alDbm^cn  +  (ty)abc^al^bmDcn]  + 

[(rz)IbcDai^bj^ck  +  (Sx)lbc^ai^bj^ck  +  (t2)£bc^ai^bjHck] 

[(rz)abcDal<5bm^cn  +  (Sz)abc^alDbB^cn  +  (tz)abc^al^bnDcn]j  • 

Performing  the  multiplication  associated  with  (6.6.9)  gives 

Afnnijk  =  S  det[Je]  |  abcWawbwc[  •  (6.6.10) 

abc 

(lxrx  +  r®ry  +  r|r|)abcDai^bj^ckDal^bB^cn  + 

(rxSx  +  lySy  +  r|s|)abcDai^bj  ^ck^alDbmfcn  + 

(rltl  +  Tyty  +  r|ti)abcDai^bj^ck^al^bmDcn  + 

(sir!  +  SyTy  +  s|r|)abc^aiDbj  ^ckDal^b*^cn  + 

(SxSx  +  SySy  +  s|s|)abc^aiDbj^ck^al^bB^cn  + 

(s|tx  +  Syty  +  sft|)abc^ai^bj^ck^alft>Bl^cn  + 

(txrx  +  t®r*  +  t|r|)abc^ai^bjDckDal^bB^cn  + 

(txSx  +  tySy  +  t|s|)abc^ai^bjDck^alDbB^cn  + 

(t|t|  +  tyty  +  t|t|)abc^ai^bjDck^al^bBDcn  ]• 

natural  conditions  exist.  Since  g  equals  n*grad$,  these  terms  cancel  along  the 
boundary  where  the  natural  conditions  apply;  the  variation  of  $  along  the 
boundaries  where  essential  boundary  conditions  applv  equals  zero.  Whence, 
these  terms  contribute  nothing  to  the  variation  of  j[$]. 
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This  expression  simplifies  after  contracting  on  the  free  indices  associated  with  the 
Kronecker  delta;  this  requires  several  steps.  First,  multiply  (6.6.10)  by  $ijk 
and  contract  on  the  indices  a,  b,  and  c  giving 

.?kAfBnijk$1jk=  (6.6.11) 

,?k  £(rfr|  +  TyTy  +  rtr!)a*nDai^Bj^nkDal  |  det[Je]  |  a«nwawBwn$?jk 
+  .E^  £(r|s|  +  rySy  +  tzSl)lbnDli^bj^nkI5b*Jdet[Je]  |  lbn^l^b^n^ijk 
+  .£  S(r|t|  +  Tyty  +  r|t|)i*cDli^»j^ckDcn  I  det[Je]  |  iBc''V]w*Wc^3jk 

ij  k  c 

+  E(s|rf  4-  sfrf  4*  sfr|)  amn^ai®  ■j^nkl^all  det[Je]  |  amnVraWBWn<j>ijk 

+  ,?k  £(SxSx  4*  SySy  +  S|sl)lbn^liHbj^nkDba|  det[Je]  |  lbnWlWbWn$ijk 

+  £  £(slt|  +  Syty  +  stt®)lnic^liDBj^ckDcni  det[Je]  |  lBCwlwiBwc$ijk 
ljk  c 

+  £  £(tlr|  +  tyXy  +  t|ri)ainri^ai^BjDnkDal|  det[Je]  |  aanwawmwn<$ijk 
ljk  a 

+  E  E(t|sl  +  tySy  4"  t|sf)lbi»£li^>jUnkUbB|det[Je]  |  lbnWlWbWn^fjk 
ljk  b 

+  £  E(t|t|  +  tyty  4-  t|t|)lBic^li5njl^ckDcn  |  det[Je]  |  lncwlwowc<|ijk- 
lj  k  C 

Using  the  properties  of  the  Kronecker  delta  while  contracting  on  the  indices  i,  j 
and  k  gives 

,?k  Afnnijk  $ijk  =  (6.6.12) 

E  (r^tj  +  Tyty  4*  r|r|)aianUaiUal j  det[Je]  |  amnwaWniWn$iBn 

i  a 

4"  £  (r|Sx  +  IySy  +  r|Sz)lbnDiiD|B  I  det[Je]  I  lbnw  1W bw n$  ibn 

lb 

+  £  (rltl  +  tyty  +  r|tf)lBcDriDcn  |  det[Je]  |  lncwlwm''*rc$iinc 

1C 

+  £  (sir!  +  SpT§  +  s|r!)aamD|jDIl|det[Je]|amnWaWBWn$|jn 

i a 

4-  £L(s|sI  4-  s§s§  4-  s|si)ibnD§jD§B|det[Je](ibnWiWbWn^fjn 
JD 

4~  £  (sftl  4"  Syty  4-  slt|)iBcDBjDcn  |  det[Je]  |  lBcWiWBWc4!jc 
J  ® 

4"  £  (tlr|  +  tyTy  4-  ttr|)aj*nl^nkUal  |  det[Je]  |  amn^a^B^n^aBk 

Jka 

4-  £  (t|sf  4-  tySy  4"  t|s|)lbnDnkU§B  |  det[Je]  |  lbnWlWbWn$5bk 
kb 

4"  £  (tit®  4"  tyty  4"  t®tz)imcDckDcn !  det[Je]  |  lBcW]WBWc^eBk- 
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The  superscripts  on  D  indicate  the  local  direction  associated  with  the 
derivative.  We  must  maintain  this  designation,  since,  in  general,  the  number  of 
nodes  in  each  of  the  three  local  directions  differ.  In  other  words,  the  quantities 
Dr,  Ds,  and  Dl  differ.  A  simplification  occurs  upon  introduction  of  the 
symmetric  array  of  coefficients, 

(Cu)amn  =  (6-6.13) 

(r|r|  +  TyTy  +  r^r^^an  I  det[Je]  I  &an^awBwn> 

(C12)5bn  H  '  6.14) 

(ris|  +  rfsy  +  r|s|)ibn  |  det[Je]  J  ibnW]WbWn, 

(C13)!»C  =  (6.6.15) 

(r|t|  +  Tyty  +  rzt|)iBc]  det[Je]  |  i»cWiw«wc, 

(C2i)aon  =  (6.6.16) 

(sir I  +  S$T§  +  sir|)anin  |  det[J  ]  |  aBn^a^B^n) 

(C22)?b„  i  (6-6.17) 

(sis!  +  S$s£  +  s|s|)lbn|det[Je]ilbnWlWbWn, 

(C23)foc  =  (6.6.18) 

(§x^x  Syty  Szt|)lnc  |  d6t[Je]  |  lmcWlWnWc, 

(C.0J..  5  (6.6.19) 

(tiri  "h  ly^y  ”1"  tzr|)amn  |  det[Je]  |  amn^a^n^nj 

(C32)ibn  =  (6.6.20) 

(tisi  +  tySy  +  t|s|)lbn  |  det[Je]  |  Ibn^l^b^tij 

and  (C33)inc  =  (66.21) 

(tit|  +  tyty  +  tztz)lac  |  det[Je]  |  IbC^IWbWc, 

into  (6.6.12),  giving 

E  Afomjk  $!jk  =  (6.6.22) 

ijk 


63 


£  (CiOimDIiDJifan  +  £  (C12)?bnDIiD|4!bn  + 

ia  ib 

£  (C13)!.cDIiD^?,c  +  £  (CI2)E,nD£jD£i4ljn  + 

ic  ja 

£  (C22) IbnDIjDg^fin  +  £  (C2 j)f BcD£j Dcn^f j c  + 
jb  jc 

£  (Cii)S.nDSkD5ifck  +  £  (C23)?bnDSkD6^!bk  + 

ka  kb 

£  (C3S)?.cD£kD^?.k. 
kc 

Reordering  the  terms  in  this  relation  yields 

£  Afanijk  $jk  =  (6.6.23) 

ijk 

£  (Cn)|anDaiDal$imn  +  ?  (Ci2)aanDalDi$ajn  + 
ia  J  a 

£  (Cl3)i»nD£lD£k<&nk  +  £  (C12)?bnDEBDii$?bn  + 

ka  ib 

£  (C22)!bnD!jD|B$!jn  +  £  (C23)!bnD|mD£k^bk  + 

j  b 

£  (Ci3)fmcDcnDTi^iBC  +  £  (C23)?mcDcnD«j$fjc  +  £  (C33)fmcDckDcn^fmk- 
ic  jc  k c 

Further  arrangements  produce 

£  Afanijk  $!jk  =  (6-6-24) 

ijk 

(Cn)lmn(D£ifc)  +  ?  (Cl2)|«n(D|^Iin)  +  £  (C13)tffl„(D‘k$£mk)]  + 

a  i  J  * 

£[£  (C12)!bn(Dl4?bn)  +  £  (C22)5bn(D|jfcn)  +  £  (C23)?bn(D£k$?bk)]Dgn  + 

b  i  J  k 

£[£  (C13)fBc(DIi^Bc)  +  £  (C23)!BC(DIj$?jc)  +  £  ( C 3 3)?BC(Dck$f mk)] D cn 4 
c  i  J  k 

Introducing 

B£|„  =  (6.6.25) 

BI5„  =  EDSjfen,  (66.26) 

and  Blau  =  EDikfck.  (6-6.27) 

k 

and  substituting  into  (6.6.24)  gives 


JThe  quantities  within  parentheses,  for  example  Dalian,  can  be  written  as  Aam 
after  performing  the  matrix  multiplication.  In  this  form  the  product  DaiAamn 
does  not  satisfy  the  rules  of  matrix  multiplication.  In  matrix  multiplication  the 
second  index  on  Dai  must  equal  the  first  index  of  Aann.  This  involves 
reorienting  the  indices  of  Dai  by  using  the  transpose.  Hence  the  product  should 

be  written  as  (Dia^Aamn,  where  the  superscript  t  refers  to  the  transpose. 
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(6.6.28) 


,?k  AWjk  #jk  - 

^(D^a)t[(Cn)^jnBa£n  +  (CnJaanBann  +  (Cl3)aMnBlan]+ 

£[(C12)?bnBI|D  +  (C22)?bnB?gI1  +  (C23)!bnB|gn]D|s+ 

^[(Cl3)?*icBInc  +  (C23)5bcB|»c  +  (C33)fBcBfgc]Dcn- 
c 

This  equation  applys  throughout  a  single  element.  Equation  (6.6.7)  summed 
over  all  elements  gives 

.?k[Af«nijk  “  ^2BfBnijk]$ijk  = 

E 

S  .?  n^nnijkFijk-  (6.6.29) 

This  matrix  expression  for  the  unknown  expansion  coefficients  $?Bn  appears  in 
the  standard  form,  Ax=B.  These  coefficients  represent  the  velocity  at  the 
nodes.  We  use  this  expression  in  a  global  node  numbering  system,  where  the 
global  nodes  lying  on  faces  shared  between  adjacent  elements  receive  a  single 
label  and  contain  a  single  value  for  each  variable,  though  several  elements  may 
contribute.  The  assembly  process,  referred  to  as  ‘direct  stiffness’,  properly 
accounts  for  the  contributions  from  several  elements. 
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APPENDIX  A 

CO-ORDINATE  TRANSFORMATION 
A.l  Global  and  Local  Co-ordinates 

We  decompose  the  given  domain  into  arbitrarily  shaped  six-sided  elements, 
we  then  transform  these  elements  from  the  global  coordinate  system  which 
applies  to  the  entire  domain  to  a  local  co-ordinate  system  within  each  element. 
The  elements  in  physical  space  transform  into  cubes  in  local  space,  with  local  co¬ 
ordinates  ranging  between  ±  1.  Since  the  physical  equations  apply  in  global 
space,  we  must  develop  a  method  for  representing  spatial  derivatives  and 
integrals  in  the  local  system.  Specifically,  we  require  a  method  to  transfer  from  a 


coordinate  system  (xi,X2,X3)  to  a  system  (£1,^2, £3). 

Assuming  an  expression  such  as 

*i  =  Xi(£h£2,6)  (A. 1.1) 

exists  between  the  two  coordinate  systems,  the  total  derivative  of  xi  appears  as 

^  =  If}*^  =  +  +  (A. 1.2) 

Expressing  this  relation  as  a  vector  by  letting  dx*  equal  the  vector  with 
components  (dxi/5£j)d£j,  changes  the  total  differentials  into 

dx1  =  +  f|^3e3>  (A.1.3) 

dx2  =  |^d£iei  +  (A. 1.4) 

and  dx3  =  ^p^d£iei  +  ^j^d£2®2  +  ^^d^e3-  (A. 1.5) 

We  require  differentials  in  vector  form  since  the  various  surface  integrals 


encountered  in  the  analysis  contain  differential  areas  pointed  in  the  surface 
normal  direction. 

A.2  Transform  Between  Global  and  Local  Coordinates 

Equation  (A.l. 2)  represents  a  valid  procedure  for  transforming  from  the 
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global  co-ordinates  (x,y,z)  to  the  local  co-ordinates  (r,s,t).  In  matrix  form  this 
transformation  equals 

fdx]  fxrXsXtlfdr]  Fdr] 

(A-2.1) 


(A. 2.2) 


dx 

Xr  XS  Xt' 

dr' 

dr’ 

dy 

= 

yr  ys  yt 

ds 

=  [J] 

ds 

dz 

Zr  Z§ 

dt 

dt 

where  the  Jacobian,  [J],  equals 

Xr  x8  Xt 

yryByt  • 

Zg  Zt_ 

Subscripts  refer  to  partial  derivatives  with  respect  to  the  local  co-ordinates.  The 


[j]  = 


inverse  transform,  whereby  the  local  co-ordinates  appear  in  terms  of  global  ones, 
equals 


dr" 

dx 

ds 

=  [J'1] 

dy 

dt 

dz 

The  inverse  of  the  Jacobian  equals 


[Ji  = 


i 

HetJJf 


ysz t-ytZs  xtz^-xgzt  xsyt-xtys 
ytzx-yrzt  xrzt-xtzr  xtyr-xryt 
yrZs-ysZr  XsZr-XrZg  Xrys-Xsyr. 


which  simplifies,  upon  introduction  of 

xj  =  yszt-ytzs,  x2  =  ytzr-yrZt,  x3  =  yrzs-yszr, 
yi  =  xtZs-xsZt,  y2  =  xrzt-xtzr,  y3  =  XsZr-xrzs, 

zi  =  xsyt-xtys,  Z2  =  xtyr-xryt,  z3  =  xrys-xsyr, 


(A.2.3) 

(A.2.4) 

(A.2.5) 


and 

to 


*  HetJJJ 

Substituting  into  the  inverse  transform  results  in 


xi  yt  z, 
x2  y2  z2 
xj  y3  z3 


dr] 

ds 

dt 


1 

astpj 


Xi  yi  Zi 
x2y2z2 
ixs  ys  z3j 


1  rdx] 
dy 


dz 


(A.2.6) 


(A. 2. 7) 


which  represents  a  convenient  method  for  computing  total  derivatives  of  local  co¬ 
ordinates  as  functions  of  global  differentials.  Partial  derivatives  of  the  local  co¬ 
ordinates  occur  when  two  of  the  global  co-ordinates  remain  constant,  or, 
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equivalently,  by  setting  one  or  more  of  the  differential  quantities  dx,  dy,  or  dz 
equal  to  zero. 

A.3  Volume  Transformation 

The  variational  form  of  the  equations  of  motion  contain  volume  integrals  in 
physical  space;  however,  the  expansion  functions  only  apply  in  transform  space. 
Consequently,  we  require  a  transformation  procedure  for  conversions  between  the 
two  systems. 

If  the  vector  dx1  represents  a  differential  length  in  the  global  system,  a 
differential  volume  equals 

dV  =  dx1*(dx3*dx3).  (A.3.1) 

Using  the  previously  derived  expressions  for  the  vectors  dx1  gives 

dV  =  cyk^d^i  f|jd£j  f|^dfk,  (A.3. 2) 

where  Cjjk  equals  the  alternating  unit  tensor.  Introducing  the  Jacobian  gives 

dV  =  det[J]dfcdfcdfc;  (A.3.3) 

where  the  determinant  of  the  Jacobian  depends  on 

det[J]d^jd^2d^3  =  €ijk^j  f|^fid£jd£k 

(AX4) 

We  require 

0  <  |  det[J]  |  <  ®  (A.3.5) 

for  a  proper  transformation.  With  (A.3.3)  we  have  a  method  for  transforming 
volume  integrals  between  the  global  and  local  co-ordinate  systems.  This 
transformation  appears  in  the  variational  statement  of  the  governing  differential 
equations,  or  whenever  a  volume  integral  in  the  global  system  appears. 

A.4  Partial  Derivatives 

The  governing  differential  equations  require  partial  derivatives  of  the 
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dependent  variables  in  terms  of  the  global  co-ordinates.  Since  the  global  co¬ 
ordinates  depend  on  the  local  co-ordinates,  and  vice  versa,  partial  derivatives 
may  appear  as 


dp  _  dpdi  ,  dpd s  ,  dpd t 

H  drdx  ^  dsdx  dtdx’ 

(A.4.1) 

dtp  _  dipdt  .  dpd s  dpd t 
dj  drdy  +  dsdj  +  dt&y’ 

(A.4.2) 

and 

dp  dpdi  ,  bpd%  ,  dpdi 

dz  wm  +  im  +  wfc- 

(A. 4.3) 

The  partial  derivatives  of  local  co-ordinates  depend  on  (A.2.7). 

the  partial  derivative  dijdx  equals 

For  example, 

dr  __  dr  _  1  _ 

cE  3x  y,z  detjjj  l' 

(A.4.4) 

Using  similar  expressions  for  the  other  derivatives  gives 

dp  xi  dp  ,  X2  dp  .  x3  dp 

W  ~  det]3fcjF  +  aei'[J]'5s  +  det  [Jf9t  ’ 

(A.4.5) 

dp  _  yi  dp  ,  72  dp  ,  y3  dp 

dy  ~  3et]Jf7r  +  JetPjtfl  +  «3?t[JJ'3t’ 

(A.4.6) 

and 

dp  Zi  dp  ,  Z2  dp  ,  z3  dp 

df  =  detppr  +  'detlJps  +  detfjJSt* 

(A.4.7) 

Expressions  such  as  these  appear  whenever  we  represent  partial  derivatives  of  the 
dependent  variables  in  transform  space. 

A.5  Partial  Derivatives  Using  Spectral  Expansions 


The  approximation  procedure  expresses  the  dependent  variables  as  a 
truncated  series  of  interpolating  functions  and  coefficients.  Therefore, 
expressions  such  as  (A.4.5),  containing  the  dependent  variable  p,  actually 
contain  a  series  of  coefficients  and  functions  which  depend  on  the  local  co¬ 
ordinates.  Our  method  uses  interpolating  functions  based  on  the  Legendre 
polynomials.  For  example,  the  one-dimensional  interpolation  function  of  order  n 
for  the  dependent  variable  p  equals 

W  (A. 5.1) 

which  contains  the  interpolating  functions,  {^k(r),  k=0,l,...,n},  and  coefficients. 
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fa.  The  derivative  of  ip  with  respect  to  the  local  variable  r  appears  as 

Dny<r)  =  J0Pkgj{&(r)].  (A.5.2) 

In  the  current  formulation,  we  compute  the  derivatives  at  the  nodes  and  not 
between  nodes.  Evaluating  the  derivative  at  ra  gives 

Dn^ra)  =  £  ^bkDaki  (A. 5. 3) 

k*0 

where  ^^k(r)]r=a=Dak  represents  the  derivative  of  the  expansion  function  fa 
with  respect  to  the  local  co-ordinate  r  at  node  ra.  Substituting  similar 
expressions  for  the  other  partial  derivatives  in  (A.4.5)  gives 

(^)abc  =  (a'et|J|)^bc  (A. 5.4) 

+  (det|j])*bc  if k^*jk^^ra^b'’ ^detfjp*1*  ifk^»jkV'i(ra)V'j(sb)Dck- 

The  superscript  e  represents  the  element  number,  while  the  subscripts  abc 

stand  for  the  spatial  node  under  consideration.  Since  the  interpolating 
polynomials  satisfy  V’i(ra)=^ai,  the  partial  derivative  expressions  reduce  to 

(G.5.5) 

(^)abc=  (jet  [Jj)^bc  f  ^bcDai  +  (^^'Jj)abc  f  ^ajcDbj  +  (det|J]^bc  f  ^abkDck, 
(|^)abc=  (g-e{|j])abc  f  ^ibcDai  +  (jgf|-j|)abc  f  ^ajcDbj  +  (ggfpj)abc  f  ^Ibk^ck, 
(^)abc=  (def[Jj)abc  f  ^lbcDai  +  (^|jj)abc  f  ^IjcDbj  +  (jg*|j])abc  f  V’abkDck- 
These  equations  represent  the  partial  derivatives  of  the  dependent  variable,  y?, 
at  the  nodal  points.  At  this  point  we  impose  no  restriction  on  the  location  of  a 

node  within  an  element,  i.e.,  the  node  may  lie  in  the  interior  or  on  the  surface  of 

an  element.  The  next  section  computes  derivatives  along  element  faces  by 
applying  (A.5.5)  to  nodes  lying  along  the  surface. 

A.6  Derivatives  Along  Element  Faces 

The  calculation  of  surface  stresses  requires  derivatives  along  element  faces. 
The  equations  for  derivatives  at  any  node  (a,b,c)  simplify  along  the  element 
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surface  since  one  of  the  indices  remains  fixed  while  the  other  two  vary.  This 
results  in  a  set  of  expressions  for  each  of  the  six  dement  faces: 

(||)abl=  (jef[j[)abl  S  £ibll>8i  +  (geqjj)abl  £  pljlDbj 

+  (getfjj)Ibl  f  (A.6.1) 

(|^)Ibi=  (aet'fjj)*bI  ?  ^*b,Dai  +  (jetfjj)*bl  ?  ^JiDbj 
+  (det  [Jp^bl  f  ^bk^lk 

(|f)!bi=  (gefpj^bl  ?  ^bl°ai  +  (cleffjp*bl  f  ^ilDbj 
+  ((fef[Xj)abl  f  V’abkDjk,  along  the  side  at  t=  - 1; 

(|f)Ijtc=  (aef  j-Jl^jtc  f  V’ljtcDai  +  (agffjpajte  ?  ^IjcDjtj 

+  (det  [  Jp*Jtc  £  ^jtkDck  (A. 6. 2) 

(f^)ajtc=  (a^Jpaitc  ?  ^jtcDai  +  (gg^Jpajtc  ?  ^IjcDjtj 

+  ^3etpJ^j‘c  k  ^tkDck 

(|f)ljtc=  (dif[Jj)ijtc  ?  ^jjtcDai  +  (g^pj)ajtc  ?  V’ajcDjtj 
+  (ggfJjj)ijtc  s  ^IjtkDck,  along  the  side  at  s=l; 

(^)tbkt=  (jet[j])^bkt  f  ^ibkt-^ai  +  (get'[T])abkt  ?  ^ajlctDbj 

+  (dgffjpabkt  |  ^IbkDktk  (A.6.3) 

(^y)abkt=  (def|Jj)abkt  ?  ^ibkt^ai  +  (det  [Jj^bkt  f  ^ajktDbj 

+  (aiftrj^bkt  f  ^bkDktk 

(^)Ibkt=  (det[Jj)abkt  f  ^ibkt^ai  +  (det  [Jpabkt  j*  ^ajktDbj 
+  (diffjj)»bkt  £  ^abkDktk,  along  the  side  t=l; 

(|f)!lc=  (difpj)alc  ?  ^ticDai  +  (ggqjpate  S  ^ajcDij 

+  (def[Jp^lc  k  ^a,k^ck  (A. 6.4) 

{^P%\c~  (deffjf)alc  f  ^ilcDai  +  (deffj])alc  ?  ^ajc^lj 

+  (det|JJ^lc  k  ^lkDck 

(|f)Ilc=  (diffjj)alc  ?  ^»!lcDai  +  (g^Jj)aic  ?  ^IjcDjj 
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(A.6.5) 


+  g^|jj)aic  s  v’likDck,  along  the  side  s=  -  1; 

(|f)1bc=  (det[J])^bc  f  ^bcDii  +  (jet  [J])^bc  ?  ^IjcDbj 

+  l  ^tbkDck 

01bc=  (ae?IJl^bc  ?  ^bcDli  +  (aetfjj)^  J  ^icDbi 

+  (detfj]^bc  f  ^bkDck 

(Ifjtbc-  (aifjjj)ibc  s  ^ibcDii  +  (a^pj)!bc  ?  ^jcDbj 
+  (gifpj)tbc  s  v>1bkDck,  along  the  side  r=  - 1;  and 

($£)!tbc=  (jlet  ["J))'tbc  f  ^ibcDiti  +  (deT[3])^tbc  j'  ^tjc^bj 

+  (<jef|jj)itbc  S  v>?tbkDCk  (A.6.6) 

(fy)  *tbc=  (aet[Jj)?tbc  f  ^bcDiti  +  (det[J]^tbc  j  ^jc^bj 
+  (detffjj^tbc  f  ^tbk‘Dck 

(ff)1tbc=  (aef  [Jj)*tbc  f  +  (det[j]^tbc  f  ^licDbj 

+  (diffjj)^bc  J  ^itbkDck,  along  the  side  r=l. 

A.7  Surface  Integrals 

Often  surface  integrals  involving  the  differential  area  vector  appear  in 


boundary  conditions  and  computations  of  surface  stresses  inter  alia.  According 
to  (A. 1.3)  through  (A. 1.5),  the  transformation  between  global  co-ordinates 
(x,y,z)  and  local  co-ordinates  (r,s,t)  equals 


dx  =  |jrdrer  + 

(A.7.1) 

dy  =  fr^01,  ^s^se®  "**  ft*^et’ 

(A. 7.2) 

and 

dz  =  -gjdrer  +  ^dsCs  +  ^dtet- 

(A. 7.3) 

In  the  local  co-ordinate  system,  the  element  surface  depends  on  two  of  the  three 
local  variables  while  the  third  remains  fixed.  For  example,  along  side  one  co¬ 
ordinates  r  and  s  vary,  while  co-ordinate  t  equals  —1.  Hence,  the 
transformation  along  side  one  appears  as 
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dx1  =  fs^ses» 

(A.7.4) 

dy1  =  |jdrer  +  Ifdseg, 

(A.7.S) 

and 

dz1  =  ^drer  +  ||dse8, 

(A.7.6) 

where  here  the  superscript  refers  to  side  one  and  not  direction  one  as  in  (A. 1.3). 


After  integrating  these  relations,  expressions  for  the  physical  co-ordinates  in 
terms  of  the  local  coordinates  along  the  surface  appear: 


xi  =  x(r,s), 

(A.7.7) 

yi  =  y(r,s), 

(A.7.8) 

and 

z1  =  z(r,s). 

(A. 7.9) 

Let  R^x^x+yJey+z^  represent  a  position  vector  describing  the  location  of 
any  point  along  surface  one.  Substituting  the  previous  expressions  for  the 
physical  co-ordinates  into  this  position  vector  gives 

R1  =  x(r,s)ex  +  y(r,s)ey  +  z(r,s)es..  (A. 7.10) 

A  differential  area  vector  in  the  normal  direction  to  the  surface  corresponding  to 
a  differential  change  in  the  local  coordinates  r  and  s  equals 

dA1  =  —  (|y*|y)drds.  (A. 7. 11) 


The  minus  sign  indicates  that  the  area  vector  points  in  the  local  co-ordinate 
direction  -  et-  We  can  write  the  two  partial  derivatives  as 


5R  5R5x„  ,  dRdy„  ,  dRdz„ 

w  ~  wtF*  +  +  • 

(A.7.12) 

and 

3R  <9Rdx_  ,  dRdy„  ,  dR dz„ 

7s  ~  dxds e*  1  '*  dzd s®** 

(A.7.13) 

Using  (A.7.10) 

reduces  these  expressions  to 

dR 

-gj  —  Xrex  +  yj€y  +  ZiOz, 

(A.7.14) 

and 

If  =  Xsex  +  ysey  +  z^. 

(A.7.15) 

Substituting  these  expressions  into  dA1  and  taking  the  cross-product  gives 
dA1  =  -  [(yrzs  -  zrys)ex  +  (zrxs  -  xrzs)ey  + 
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(xrys  -  yrxs)ez]drds, 


(A.7.16) 


or,  using  the  simplified  notation  of  (A.2.5) 

dA1  =  -  (x3ex  4-  y3ey  +  z3ez)drds.  (A. 7.17) 

Similar  expressions  for  sides  two  through  six  equal 

dAJ  =  (xiex  +  y2®y  +  zjez)drdt,  (A.7.18) 

dA3  =  (x3ex  +  y^  +  z^drds,  (A.  7. 19) 

dA4  =  -<X2ex  +  y2ey  +  zje^drdt,  (A.7.20) 

dA5  =  — (xiex  +  yiey  +  ZiCz)dsdt,  (A. 7.21) 

and  dA®  =  (xje*  +  yiey  +  ziez)dsdt.  (A.  7.22) 


Expressions  such  as  these  will  appear  whenever  we  compute  surface  integrals. 
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APPENDIX  B 


SURFACE  FORCE  AND  MOMENT  CALCULATION 
B.l  Background 

Fluid  flowing  past  solid  boundaries  develops  surface  stresses  which  lead  to 
surface  forces  and  moments  after  integration  over  the  surface  area.  The  elements 
located  adjacent  to  these  surfaces  contain  flow-field  information  necessary  for 
surface  stress  calculations  Let  dA  equal  a  differential  area  on  the  element 
surface  with  unit  normal  vector  n.  Let  the  force  per  unit  area  due  to  the  fluid 
stresses  on  this  surface  equal  t(n).  Then,  the  total  force  on  a  finite  surface,  dU, 
results  from  integrating  the  contributions  from  all  the  differential  elements: 

F=  ff  t(n)dA.  (B.1.1) 

Similarly,  the  total  moment  about  a  given  point  due  to  the  surface  forces  equals 

M  =  ff  x»t(n)dA,  (B.1.2) 

dQ 

where  x  represents  the  distance  between  the  surface  element  and  the  point 
about  which  the  moment  is  taken.  Let  <7ji  denote  the  ith  component  of  t(j), 
and  let  t(n)j  denote  the  ith  component  of  t(n).  Whence 

t(n)i  a  crjinj,  (B.l. 3) 


where  the  stress  tensor  equals 

0  11  <712  cr  13 

£  =  (721  022  023  •  (B.1.4) 

(731  ^32  033. 

Since  t(n)  represents  a  vector  and  n  represents  a  unit  vector  independent  of 


the  <7ji,  we  require 


t(n)  =  n*£. 


(B.l. 5) 


This  equation  represents  the  force  per  unit  area  due  to  fluid  stresses  on  a  surface 
with  unit  normal  vector  n.  When  this  expression  replaces  t(n)  in  the  force  and 


moment  integrals,  terms  such  as  n*cndA  and  x*(n-ff)dA  arise.  Since  dA 
represents  a  scalar,  we  may  write  these  expressions  as  ndA-£  and  x*(ndA-<r). 
We  developed  expressions  for  ndA=dA  along  each  of  the  six  element  faces  in 
appendix  A.  We  still  need,  however,  an  expression  for  the  stress  tensor. 

B.2  Incompressible  Stress  Tensor 


In  Cartesian  co-ordinates  the  diagonal  components  of  the  incompressible 
flow  stress  tensor  equal 


n  t\ 

0-XX  =  -  P  +  <7yy  =  -  p  +  2^-, 

and 

0ZZ  =  —  P  +  2fJgp 

(B.2.1) 

while  the  off  diagonal  terms  equal 

fd u  ,  3v> 

Gxy  -  <Tyx-  /A^y  +  ^)> 

(B.2.2) 

<Txz  =  <?zx  =  /<f|  +  ff), 

(B.2.3) 

and 

<7yZ  =  (Jzy  =  /i(|j  +  ff  )• 

(B.2.4) 

B.3  Surface  Force 

As  indicated  in 

§B.l,  the  force  over  a  surface  30  due  to  the  fluid  stresses 

equals 

F  =  ff  t(n)dA  =  ff  n-odA. 

*"30  *"30  ~ 

(B-3.1) 

The  components  of  the  force  in  the  global  (x,y,z)  system  equal 

Fx  —  ff  (oxx^x  +  CyxUy  +  0zxflz)dA, 

•"30 

(B.3. 2) 

Fy  —  (fl^y^X  4"  (TyyHy  *4"  ffj5ynz)dA, 

(B.3. 3) 

and 

Fz  =  (axZnx  4-  Cyzny  +  (Tzznz)dA. 

30 

(B.3.4) 

Since  nxdA,  nydA, 

and  nzdA  represent  the  projection  of  dA 

onto  the  planes 

perpendicular  to  the  global  co-ordinate  axes,  we  may  write  nxdA=ex*dA, 
nydA=ey*dA  and  nzdA=ez-dA.  Substituting  these  expressions  into  the  force 
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components  gives 


Fx  =  ff  (<TXXeX  +  (Tyjfiy  + 

(B.3.5) 

Fy  =  //  (ffxyex  +  Oyy€y  +  <7zy€z)*dA, 

',dCl 

(B.3.6) 

and 

Fz  =  ff  (<Tt2.£x  +  0"yz6y  +  (7zZez)*dA. 

JJdQ 

(B.3.7) 

The  forces  along  side  one,  where  dAJ=  —  (xsex+yjey+zsezjdrds,  equal 

Fxl  =  -  ff  (<rxxx3  +  owy3  +  trzxz3)drds. 

(B.3.8) 

Fy1  =  —  ff  ( G xyx 3  +  OyySl  'F  &zyZz)dids, 

JJdW 

(B.3.9) 

and 

Fz1  =  -  ff  (axzx3  +  ayzy3  +  Czz z3)drds. 

JJ8U  i 

(B.3.10) 

On  the  local  level,  the  surface  of  integration  reduces  to  a  square  with  sides 
ranging  between  ±1.  Using  the  spectral  collocation  approximation  allows  an 
expression  for  the  force  Fxl  as 

Fxie=  -£  /  lf\<rxxX3  +  <7yxy3  + 

^zx23)iji^i(r)^j(s)dxds,  (B.3.11) 

where  the  expression  indicates  that  the  terms  within  the  parentheses 

depend  on  the  local  node  ijt  of  element  e.  Integrating  over  the  element  surface 
gives 

Fxie  =  -— £  (&xxx3  +  &yxj3  +  ^zxZ3)ijiWiWj.  (B.3.12) 

In  a  similar  manner,  the  y  and  z  components  appear  as 

Fyie=  “£  ( O' XyX 3  +  CTyyy 3  +  CTzyZ3)ljlWiWj  (B.3.13) 

and  Fzie=  — E_  (.txzx3  +  nyzys  +  OzzZs^jiWiWj.  (B.3.14) 

The  forces  on  faces  two  through  six  equal 

Fx2e  =  E^  (crxxX2  +  ^yxj2  +  0zxz2)ijtkWiWk, 

Fyje  =  E^  (<7XyX2  +  0yyy2  +  <7zyz2)ijtkwiwk> 

and  =  E  ((7X zx2  +  0yzy2  +  Ozz^ljtkWiWk,  (B.3.15) 

ik 
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on  face  two,  where  dA^x^x+y^y+z^jdrdt; 

Fx3e  =  E  (^xxX3  +  ffyxy3  +  aZxZ3)?jktWiWj, 

Fy3C  =  ?j  (JxyXs  +  ayyyi  +  ®*j*»)tiktWiWj, 
and  Fz3e  =  E  (<rxzx 3  +  tfyzya  +  trZaZ3)?jktWiWj,  (B.3.16) 

on  face  three,  where  dA=(x3ex+y3ey+Z3ez)drds; 

Fx4®  =  E  (<Ttx*2  +  0yjj2  +  ffzx*2)ilkWiWk, 

Fy4®  =  Sk  {axyX2  +  ayyy2  +  trtyz2)!ikWiWk, 
and  Fz4e  =  (<rxzx2  +  <ryz y2  +  crzzz2)tikWiWk,  (B.3.17) 

on  face  four,  where  dA=  -  (x2ex+y2ey+z2ez)drdt; 

Fx5e  =  ?k  (<7xxxj  +  (Tyjji  +  crzxZi)tjkWjWk, 

Fy5C  =  ?k  (^xyX!  +  ayyyi  +  ffzyZOtjkWjWk, 

and  Fzse  =  (^axt  +  ayzyj  +  azzzOljkWjWk,  (B.3.18) 

on  face  five,  where  dA=  —  (x  iex+y  iey+z  iez)dsdt ;  and 

Fx«e  =  ?k  (axxXi  +  cryx 7 1  +  azxzOftjkWjWk, 

Fy86  =  Ek  (ffxyXl  +  (Tyyjl  +  <7jyZl)?tjkWj Wk, 

and  Fz«e  =  E^  (axzXj  +  <ryzyi  +  <rzzzi)?tjkWjWk,  (B.3.19) 

on  face  six,  where  dA=(xiex+yiey+Ziez)dsdt.  Note  that  when  any  of  the 
subscripts  ijk  equals  1,  we  refer  to  the  node  at  the  face  corresponding  to  rst 
equal  to  —1,  and  when  any  of  the  ijk  equals  it,  jt,  or  kt,  we  refer  to  the 
node  at  the  face  corresponding  to  rst  equal  to  1.  Also,  the  basis  functions 
integrated  over  the  local  co-ordinate  domain,  Wi=J*  j^i(r)dr,  Wj=J*  ^j(s)ds, 
and  Wk=J"  |V<i(t)dt,  differ  when  the  number  of  nodes  in  the  r,  s,  or  t 
directions  differ;  consequently,  we  must  remember  that  wi  refers  to  the  r- 
direction,  wj  refers  to  the  s-direction,  and  Wk  refers  to  the  t-direction. 

B.4  Surface  Moment 
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The  moment  along  a  surface  due  to  the  fluid  stresses  appears  as 

M  —  ff  x*n*ndA,  (B.4.1) 

wHre  the  moment  "arm"  equals  x=Xx^x+Xyey+Xzfiz-  Substituting  into  (B.4.1) 
and  performing  the  cross  product  yields  the  moment  components  in  the  global  co¬ 


ordinates  as  (B  4  2) 

Ml  =  ff  [Xyi^xz^x  +  (Tyz^y  +  (Tzzfiz)  Xz(^xy^x  +  CyyXiy  +  <Tzynz)]dA, 

"dn 

My  =  ff  [Xzi^xx^x  4*  Oyx&y  +  OzxBz)  ~  Xx(^xzAx  +  (Tyz&y  4*  Ozznz)JdA, 

en 

Mz  —  ff  [XxC^Xynx  +  CTyylly  +  Czyllz)  Xy(^XXnx  4*  OyxP-y  4*  OzxnZ)]dA. 
Using  the  same  notation  as  in  the  surface  force  calculation  gives 

Mx  =  JT  [Xy(<7xzex  +  (7yzey  4"  O'zzBz)  —  Xx(^xy®x  d"  PyyCy  +  &Z yCz)]  *  d  A, 

'"dn 

My  =  ff  [Xz{^xx^x  4-  f'yxfiy  +  PzxPz)  ~  Xx(&X?£z  +  (TyzOy  +  <Tzz&r)] -dA, 

•'•'an 

M Z-  ff  [Xx(ffxye*  +  ayyZy  +  Oz y*z)  -  Xy(^xzex  4-  ayxey  +  azxez)j-dA. 

(B.4.3) 


For  example,  the  x-component  evaluated  along  side  one,  where 
dA=— [xsex-fyaey+zaejdrds,  equals 

Mxie  =  -/_j/_j[X y(^xzX3  4-  cryzy3  +  <rzzz3)  -  Xz(^xy*3  +  OyyYz  +  <rzyz3)]edrds. 
Introducing  the  spectral  expansion  functions  gives 

Mxl®=  — E_  J  |[Xy(^xzX3  4-  cryzy3  (B.4.4) 

+  Oi&i)  -  Xz{°x yX3  +  ayyy3  +  <rzyz3)]ijii/>i(r)^j(s)drds. 

After  integration  we  obtain 

Mxie-  -E.  [Xy(^xzX3  4-  tryzy3  +  <rzzz3) 

-  Xz(^xyX3  +  <ryyy3  +  <rzyz3)]tjiWiWj.  (B.4.5) 

In  a  similar  manner,  the  y  and  z  components  appear  as 
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Myie  =  — £.  [Xz(<7xxx3  +  PyxY*  +  0"zxZ3)  ~  Xx(ffxsXs  +  ^yzys  +  tfzaZS^ijiWi^j} 

Mzie  =  — £.  [Xx(^xyx3  4*  OyyYz  +  ffZyZ3)  —  Xy(°xxx3  4*  OyxJS  +  OzxZ3)]ijlWiWj. 

The  moments  along  faces  two  through  six  equal 

Mxje  —  [Xy (^xz^  +  <ryzy2  +  OztZ?)  —  Xz(<7xyx2  +  OyfSl  +  <7zyZ2)]  tjtkWjWk, 

Myje  =  £  [Xz(^xxx2  +  ayzY  2  +  0zxZ2)  —  Xx(°’xzx2  +  <7yzY2  +  0zzZa)]ijtkwiwk> 

ik 

M/  =  £  [Xx(^ryx2  +  Gyyji  +  0zyZ2)  —  Xy(^xxx2  +  Oyxj2  +  <7zxZ2)]ijtkwiwk, 
ik 

along  face  two,  where  dA=(x2ex+y2ey4-z2ez)drdt;  (B.4.6) 

Mxje  =  £  [Xy(°‘xzx3  4-  GyzYz  4-  azzzl)  ~  Xz(<7xyx3  +  ffyy7z  +  0zyz3)]ijktWiWj, 

My3e  =  £  [Xz(^xxx3  +  PyxYl  +  &zx%z)  ~~  +  Pyzjl  +  ^zzZs)]  ijktw iwj) 

Mz3e  =  £  [Xx(<7xyx3  +  O’yyYl  +  <XZyZ3)  —  Xyt^xx^a  +  CTyxJZ  +  OzxZs)] fjktw iw j , 

along  face  three,  where  dA=(x3ex+y3ey4-Z3ez)drds;  (B.4.7) 

Mx4e  =  — £^  (Xy(ffxzX2  +  <TyzY2  +  <7^2)  —  Xz(°xyx2  4-  0yyj2  4-  0zyZ2)]ilkwiwk, 

My<e  =  -£  [Xz(<7xxX2  4-  ffyxy2  +  OzXZ2)  -  Xx(<7xzX2  4-  0ryZy2  4-  <7zzZ2)]?lkWiWk, 

ik 

Mz<e  =  — £^  [Xx(<XxyX2  4-  Cyyy2  4-  CTzyZ2)  —  Xy(axxx2  4-  Pyx? 2  +  <7zxZ2)]  llkWiWk, 

along  face  four,  where  dA=  —  (x2ex4-y2ey4-z2ez)drdt;  (B.4.8) 

Mxse  =  — £^  [Xy(^xzxl  4-  GyzY I  4-  ffzzZ  1)  —  Xz(axyxl  4-  PyyYl  4-  0’2yZi)]fjkWjWk, 
My5e  =  — £^  [Xz(<7xxxi  4-  <ryxy  1  4-  cr zxz  1)  —  XxC^xzXi  4-  GyzY  1  4-  ffzzZi)]?jkwjwk> 
Mz5e  =  -?k  [x*(<7XyXi  4-  oyyyi  4-  azyzi)  -  XyC^xxXi  4-  GyxYi  4-  <rzxzi)]?jkWjWk, 
along  face  five,  where  dA=  —  (xiex+yiey+ziez)dsdt;  (B.4.9) 

Mx®e  =  £k  [Xy(<7xzXl  4-  OyiYx  +  O’zzZi)  -  Xz(^xyxl  4-  (TyyYi  4-  <7ZyZ  1)]  it  j  kw  j  wk , 
Myfle  =  £^  [Xz(<7Xxxl  4-  PyxYl  4-  ffzxZi)  —  XxC^xzXl  4-  <7yzyi  4-  0zzZi)]itjkwjwk, 
Mz®e  =  £k  [Xx(<7xyXi  4-  <7yyyi  +  CTzyZi)  —  Xy(CTxxxl  4-  ^yxy  1  4-  <7zxZi)]itjk'*rj"wk, 
along  face  six,  where  d A=  (x iex4-y  iey+z  ^dsdt .  (B.4.10) 
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APPENDIX  C 


IMPERMEABLE  WALL  BOUNDARY  CONDITIONS 
0.1  Physical  Basis 

Along  an  impermeable  wall  the  normal  fluid  velocity  must  equal  the  normal 
component  of  the  wall  velocity,  while  the  tangential  velocity  may  vary.  We 
express  this  as 

n-V  =  n-Vwaii-  (€.1.1) 

Another  condition  on  impermeable  wall  boundaries  requires  zero  tangential 
stress,  or,  equivalently,  that  the  stress  points  normal  to  the  surface.  In  other 
words,  the  cross  product  of  the  surface  unit  normal  and  the  local  stress  vector 
equals  zero.  Thus,  an  impermeable  wall,  unlike  a  solid  wall,  cannot  support  a 
tangential  stress.  In  symbolic  form  this  boundary  condition  equals 

n*(a*n)  =  0,  (C.1.2) 

where  n  represents  the  surface  outward  unit  normal  vector,  and  a  represents 
the  total  stress  tensor  consisting  of  the  sum  of  the  isotropic  and  deviatoric  terms. 
We  must  prove  that  the  combination  of  these  two  boundary  conditions  equals  the 
single  condition  given  by 

(n-V)V  =  0.  (C.2.1) 

The  derivation  of  this  relation  follows. 

C.2  Mathematical  Derivation 


Begin  by  expressing  the  unit  normal  as 


and  the  velocity  as 


n  —  UxCx  ~r  Uy€y  4"  nzCfc 

V  =  Uex  +  Vey  +  We*. 


The  inner  product  n>V  gives 


(C.2. 2) 

(C.2.3) 
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n-V  =  Unx  4-  Vny  +  Wnz. 


(C.2.4) 


Setting  this  equal  to  the  normal  velocity  of  the  boundary  gives 

Unx  4-  Vny  4-  Wnz  =  n-Vwall-  (C.2.5) 

The  stress  vector  on  a  surface  with  normal,  n,  equals  t(n)  which  appears  as 

t(n)  =  ff-n  =  (C.2.6) 

(^XXnX  4*  (TyxJXy  +  CzX^z)^X  + 

(cxynx  ■+•  Oyy ny  +  azynz)ey  + 

(crxznx  4-  cryzny  +  <7’zzUz)ez 

according  to  appendix  B.  The  cross  product  of  this  relation  with  the  unit  normal 
equals 

nxt(n)  =  n*£*n  =  (C.2.7) 

[uy(cxznx  +  cryzny  +  ffzz^-z)  —  ^z(  ^iyUx  +  0yyUy  +  azynz)]ex  + 

[nz(<7ixni  +  <TyxUy  4-  crzxnz)  —  nx(crxznx  4-  cryzUy  4-  ffzZnz)]ey  4- 
[nx(0xyUx  +  (TyyUy  4*  <7zynz)  —  Ey(axxnx  +  OyxUy  +  ffzxUz)]^- 
Using  (B.2.1)  through  (B.2.4)  for  the  stress  components  gives  the  following 


relations  for  the  three  components  of  the  above  expression: 

(n*£-n)x  =  !  +  A® yny(fy  4  f^)  +  ~  P)  ~ 

^xUz(|y+|—)  -  imynz(2jj  -  p)  -  /*nznz(|y4-|y),  (C.2.8) 

(n*£-n)y  =  /mxnz(2|y~p)  4-  /myn z(|yfJ^)  4-  — 

Mnxnx(|^+|^)  -  MnyHx(|y+|y)  -  -  p),  (C.2.9) 

(n*£-n)z  =  /mxnx(^y.4 )  4-  /mynx(2^- —  p)  4-  MQzBx(| y  +  )  — 

Amxny(2|^  -  p)  -  /xnyny(|^4-|y)  -  /rnzny(^4-^).  (C.2.10) 

The  pressure  terms  cancel  in  each  of  these  expressions  giving 


(n*£*n)x  -  /mxny(|Y  +  f^f)  +  /xnyI1y(|y-  +  |y)  +  ^ny(2lir  ~  2fy)  ~ 

i«nxnz(|y  +  ^)  - /mznz(|j  4- |y-),  (C  2.ll) 


8.1 


(nx£*n)y  =  imxnz(2j~  -  2jj-)  +  /mynz(|y  +  |^)  +  +  |j)  - 

+  It)  “  ^yn*(|y  +  fj).  (C.2.12) 

(n*£-n)z  =  /mxnx(|y  +  |j)  +  ^nynx{2|y  -  2^)  +  /mznx(|y  +  |^)  - 

W>.v(|t  +  If)  ~  W>4f  +  Jj)-  (C.2.13) 

After  further  rearrangement  we  obtain 

(n*£*n)x  =  -fmz(  »x|y  +  ny|j  +  nz|y)  +  +  ny|^  +  nz|^— ) 

+  ^flxlx  +  ny|y  +  nz^)(ny^  ”  n*V),  (C.2.14) 

(n*£-n)y  =  — ^nx(nx|y  +  ny|j  +  +  /mz(nx|^'  +  nyf^  + 

+ /x(r^+ ny^-  + nas^)(nzU-nxW),  (C.2.15) 

(nx£-n)z  =  -/my(nx||  +  n^  +  nz|^)  +  ^nx(nx|“  +  n^  +  n^) 

+  faxfe  +  % |y  +  nz^)(nxV  -  nyU).  (C.2.16) 

The  partial  derivatives  of  the  unit  normal  vector  with  respect  to  the  spatial  co¬ 
ordinates  differ  from  zero  on  a  curved  surface.  Hence,  we  cannot  interchange  the 
derivative  operator  in  the  first  two  terms  of  each  expression  with  the  components 
of  the  normal  vector  (e.g.,  *  ^(nxU  +  nyV  +  nzW)). 

However,  the  derivative  may  move  outside  the  components  of  the  normal  vector 
if  we  only  consider  surfaces  with  zero  curvature.  When  we  restrict  ourselves  to 
this  case,  the  above  relations  may  appear  as 

a  q 

(n*a-n)x  =  -/m^nxU  +  nyV  +  nzW)  +  /my|-(nxU  +  nyV  +  nzW) 

+  ^n*Jx  +  Dyly  +  nz!z )(nyW-nzV),  (C.2.17) 

<9  a 

(nx£-n)y  =  -;mr^(nxU  +  nyV  +  nzW)  +  /m^njU  +  nyV  +  nzW) 

+  +  &y|^  4-  nz|^)(nzU  -  nxW),  (C.2.18) 

a  a 

(nx£*n)z  =  ~^ny^(nxU  +  nyV  +  nzW)  +  nxU  +  nyV  +  nzW) 

+  +  Uylfy  +  n^)(n*v  _  ny^)-  (C.2.19) 

a  a  a 

Introducing  n*  V=nxU+nyV+nzW  and  n«  V^n^ t  into  these 
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expressions  yields 

(nx£-n)x=  +  Mtt*V(nyW-ii.V),  (C.2.20) 

(n*£-n)y  =  -^(nxl^-n^n-V  +  pa.V(nzU-nxW),  (C.2.21) 

(nx£-n)z  =  -fj{nyfe-nJI}y)n-V  +  fm-V(nxV-nyU).  (C.2.22) 

After  introducing  the  cross  product,  the  three  components  of  (n*a-n)  reduce  to 
a  single  vector  expression  given  by 

(n*£-n)  =  ^((nxV)n-V  +  (n*V)n«V].  (C.2.23) 

This  equation  also  appears  as 

(n*£-n)  =  /m*(V(n-V)  4-  (n-V)V].  (C.2.24) 

Introducing  n-V=n-Vwaii  gives 

(n«£-n)  =  /m*[V(n-Vwaii)  4-  (n*V)V|.  (C.2.25) 


For  a  general  curved  boundary  moving  with  velocity  Vwan,  the  quantity 
V(n- Vwaii),  does  not  equal  zero.  However,  since  we  are  limiting  our 
consideration  to  flat  boundaries,  this  term  does  equal  zero.  In  this  case  the 
boundary  condition  reduces  to 

(n*£-n)  =  /m*[(n-V)V].  (C.2.26) 

Since  the  original  boundary  condition,  n*£*n,  equals  zero,  we  require 

/m«[(n*V)V]  =  /*(n-V)(n*V)  =  0.  (C.2.27) 

This  equation  is  satisfied  when  the  unit  normal  lies  parallel  to  the  velocity  vector 
giving  n*V=0.  This  seems  like  a  most  unusual  condition  and  we  assume  that  it 
does  not  hold  for  most  boundaries.  The  other  possibility  requires 

(n*V)V  =  0,  (C.2.28) 

and  represents  the  desired  form  of  the  boundary  conditions  given  in  §C.l.  This 
proves  that  the  two  impermeable  wall  boundary  conditions  reduce  to 
homogeneous  natural  conditions  when  we  enforce  (1)  zero  surface  curvature  to 
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ensure  V(n*  VWaii)=0,  and  (2)  a  surface  with  unit  normal  not  parallel  to  the 
local  velocity. 
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